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INTRODUCTION: 


A  tracked  free-hand  ultrasound  system  is  ideal  for  guiding  many  radiotherapy 
procedures,  as  it  can  be  performed  in  the  treatment  room.  In  partial  breast  irradiation, 
the  lumpectomy  cavity  should  be  localized  during  the  treatment  course.  While  many 
structures  in  the  breast  look  similar  to  the  lumpectomy  cavity  in  ultrasonography,  the 
scar  tissue  around  the  cavity  is  hard  and  can  be  visualized  by  ultrasound 
elastography.  To  be  clinically  successful,  the  elastography  method  needs  to  be  robust 
to  the  sources  of  decorrelation  between  ultrasound  images,  specifically  fluid  motions 
inside  the  cavity,  change  of  the  appearance  of  speckles  caused  by  compression  and 
physiologic  motions,  and  out-of-plane  motion  of  the  probe.  We  presented  a  novel 
elastography  technique  that  was  based  on  Dynamic  Programming  (DP)  and  analytic 
minimization  of  a  regularized  cost  function.  The  cost  function  incorporates  similarity 
of  RF  data  intensity  and  displacement  continuity,  making  the  method  robust  to 
decorrelation  noise  present  throughout  the  image.  We  exploited  techniques  from 
robust  statistics  to  make  the  method  resistant  to  large  decorrelations  caused  by 
sources  such  as  fluid  motion.  The  analytic  displacement  estimation  worked  in  real¬ 
time.  We  further  used  the  tracked  data,  used  for  targeting  the  irradiation,  for 
discarding  frames  with  excessive  out-of-plane  motion.  In  this  method,  we  introduced 
an  “optimal  frame  selection  function”  which  selected  the  best  frames  -using  the 
tracking  data-  to  generate  high  quality  elasticity  images.  In  addition,  we  have 
extended  the  ID  method  to  2D,  meaning  that  we  can  calculate  tissue  displacements 
more  completely.  The  2D  method  follows  the  path  of  the  ID  method  in  optimizing 
regularized  cost  functions  and  is  also  real-time.  And  finally,  we  have  developed  a  3D 
elasticity  imaging  technique  based  on  DP  and  analytic  minimization  of  a  cost 
function  for  volumetric  imaging  of  the  elasticity. 

We  also  introduced  Kalman  filtering  for  calculating  strain  images.  The  Kalman  filter 
allows  estimating  low  variance  strain  images,  while  it  doesn’t  eliminate  boundaries 
by  over-smoothing  the  strain  images.  A  major  contribution  has  also  been  made  in 
utilizing  multiple  ultrasound  images  to  calculate  the  elasticity  image.  Displacement 
estimation  is  an  essential  step  for  ultrasound  elastography  and  numerous  techniques 
have  been  proposed  to  improve  its  quality  using  two  frames  of  ultrasound  RF  data. 
We  introduced  a  technique  for  calculating  a  displacement  field  from  three  (or 
multiple)  frames  of  ultrasound  RF  data.  To  calculate  a  displacement  field  using  three 
images,  we  first  derive  constraints  on  variations  of  the  displacement  field  with  time 
using  mechanics  of  materials.  These  constraints  are  then  used  to  generate  a 
regularized  cost  function  that  incorporates  amplitude  similarity  of  three  ultrasound 
images  and  displacement  continuity.  We  optimize  the  cost  function  in  an  expectation 
maximization  (EM)  framework.  Iteratively  reweighted  least  squares  (IRLS)  is  used 
to  minimize  the  effect  of  outliers.  We  show  that,  compared  to  using  two  images,  the 
new  algorithm  reduces  the  noise  and  eliminates  ambiguities  in  displacement 
estimation.  The  displacement  field  is  used  to  generate  strain  images  for  quasi-static 
elastography. 
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BODY: 


We  are  very  excited  to  report  the  novel  achievements  of  this  research  effort.  The 
detailed  Statement  of  Work  tasks  and  a  description  below  each  task  is  provided  next. 
The  SOW  tasks  are  in  blue.  A  description  of  the  research  effort  follows  each  SOW 
item  in  black. 


1 .  Obtain  ultrasound  (US),  US  elastography  (USE)  and  CT  scans  before  the 

start  of  the  radiotherapy. 

a.  Develop  a  tracked  US  system  for  data  collection  (month  1).  We  will 

use  Polaris  optical  tracker,  already  available  in  the  ERC-CISST  for  tracking.  Polaris 
markers  will  be  attached  to  the  US  probe  and  the  US  probe  will  be  calibrated  with 
respect  to  the  trackers. 

Polaris  optical  tracker  gives  very  accurate  displacement  measurements.  However, 
they  require  line  of  sight  (as  they  are  optical)  and  also  they  take  longer  time  to  set  up 
in  a  CT  imaging  room.  Because  of  the  time  constraints  (we  had  a  short  amount  of 
time  to  set  up  the  tracking  device  in  the  CT  room  before  the  patient  shows  up),  we 
instead  used  magnetic  trackers.  Magnetic  trackers  are  not  as  accurate  as  Polaris  and 
suffer  from  distortion  of  magnetic  field  due  to  presence  of  metals.  However,  they  do 
not  require  line  of  sight  and  also  easier  to  set  up  and  also  are  inexpensive.  We  started 
by  using  a  Polaris  camera  to  collect  position  data.  However,  setting  up  the  camera 
was  cumbersome  and  the  limited  time  that  we  had  in  the  CT  room  did  not  allow  the 
required  preparations.  Therefore,  we  opted  for  a  magnetic  tracker  and  continued  data 
collection  using  a  magnetic  tracker. 


b.  Collect  US  data  from  patient  before  the  PBI  treatment  at  the  same 

time  that  CT  is  collected  (months  2-14). 

We  have  collected  data  from  16  patients  so  far.  All  data  was  acquired  in  the  CT 
room  while  the  patient  was  lying  down  on  the  CT  bed.  The  patient  did  not  move 
after  ultrasound  data  was  acquired  until  her  CT  scan  was  finished.  We  have  acquired 
both  3D  ultrasound  data  (freehand  data,  no  3D  probe  is  used)  and  palpation  data  (for 
elastography).  More  details  of  the  system  are  in  the  attached  paper  (MICCAI 
conference,  34%  acceptance  rate,  MedLine  and  ISI  indexed  and  listed).  I  received 
travel  award  from  MICCAI  (only  7  students  received  this  award)  to  present  this 
work. 

Our  approach  introduces  minimal  divergence  from  the  original  workflow  of  PBI 
treatment.  We  have  an  approved  institutional  review  board  (IRB)  protocol  to  obtain 
B-mode  and  strain  images  from  patients  who  undergo  lumpectomy.  We  acquire  the 
ultrasound  data  when  the  patients  return  for  the  CT  scan  four  weeks  after  the 
surgery.  The  data  includes  tracked  B-mode  images  scanned  over  the  lumpectomy 
bed,  tracked  real-time  strain  images,  and  RF  dara  synchronized  with  the  tracking 
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information  for  off-line  processing.  We  have  devised  a  data  collection  system  for 
this  purpose  shown  in  Figure  1 . 

Before  acquiring  the  CT  scan,  four  CT-compatible  fiducials  are  placed  around  the 
scar  of  the  surgery  on  patient’s  breast.  These  fiducials  are  additional  to  the  ones  that 
are  commonly  placed  on  patient’s  chest  and  the  resting  foam  bed.  The  foam  bed, 
shown  in  Figure  1  maintains  the  configuration  of  patient’s  body  during  the 
treatment.  The  extra  fiducials  are  used  to  locally  align  the  CT  and  ultrasound  data, 
and  the  regular  ones  are  used  for  targeting  at  the  time  of  irradiation. 


Polaris  Camera 


philips 


Ultrasound  Machine 


Data  fusion  and  display 


Patient  rests  here 


L _ T 


Tracked  transducer 


Figure  1.  The  Polaris  camera  (NDI,  Waterloo,  Canada)  tracks  the  3D  orientation  of 
the  ultrasound  probe  by  detecting  the  4  bright  spheres  attached  to  the  ultrasound 
transducer. 


The  patient  remains  in  the  same  resting  position  after  the  CT  scan.  As  the  first  step, 
the  sonographer  digitizes  the  location  of  the  fiducials  with  a  calibrated  stylus  by 
simply  touching  each  of  the  four  fiducials  with  the  stylus  in  a  certain  order 
prescribed  in  our  protocol.  Afterwards,  the  sonographer  sweeps  the  ultrasound 
transducer  over  the  lumpectomy  bed  collecting  at  least  500  B-mode  images.  Then 
within  a  few  seconds,  the  program  constructs  the  ultrasound  volume  from  the 
collected  images.  Next,  the  sonograhper  obtains  individual  tracked  strain  images 
from  the  areas  of  interest  by  gently  moving  the  ultrasound  transducer  up  and  down. 
The  visualization  program  demonstrates  the  flying  strain  image  over  the  ultrasound 
volume  (Figure  2).  Depending  on  each  case,  we  acquire  a  minimum  of  200  strain 
images.  Lastly,  we  record  sequences  of  RF  signals  as  the  sonographer  compresses 
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and  decompresses  the  tissue.  The  RF  data  synchronized  with  the  tracking 
information  facilitates  more  elaborate  strain  imaging  techniques. 


Figure  2.  The  visualization  software  shows  3D  ultrasound  and  strain  images. 


To  generate  3D  ultrasound,  we  have  developed  a  novel  method  that  uses  the 
ultrasound  scatterers  to  create  a  high  quality  3D  volume.  Out-of-plane  motion  in 
freehand  3D  ultrasound  can  be  estimated  using  the  correlation  of  corresponding 
patches,  leading  to  sensorless  freehand  3D  ultrasound  systems.  The  correlation 
between  two  images  is  related  to  their  distance  by  calibrating  the  ultrasound  probe: 
the  probe  is  moved  with  an  accurate  stage  (or  with  a  robot  in  this  work)  and  images 
of  a  phantom  are  collected,  such  that  the  position  of  each  image  is  known.  Since 
parts  of  the  calibration  curve  with  higher  derivative  gives  lower  displacement 
estimation  error,  previous  work  limits  displacement  estimation  to  parts  with 
maximum  derivative.  In  this  paper,  we  first  propose  a  novel  method  for  exploiting 
the  entire  calibration  curve  by  using  a  maximum  likelihood  estimator  (MLE).  We 
then  propose  for  the  first  time  using  constrains  inside  the  image  to  enhance  the 
accuracy  of  out-of-plane  motion  estimation.  We  specifically  use  continuity 
constraint  of  a  needle  to  reduce  the  variance  of  the  estimated  out-of-plane  motion. 
Simulation  and  real  tissue  experimental  results  are  presented.  The  published  paper 
attached  at  the  end  of  the  report  titled  “Novel  Reconstruction  and  Feature 
Exploitation  Techniques  for  Sensorless  Freehand  3D  Ultrasound”  thoroughly 
explains  this  technique. 
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c. 


Register  US  to  the  CT  (months  2-14). 


We  advanced  significantly  and  proposed  exciting  new  avenues  for  further  work  in 
developing  novel  techniques  for  elastography,  which  is  the  major  contribution  of  this 
research  proposal.  Therefore,  we  have  not  had  enough  time  to  focus  on  this  task. 
However,  we  have  performed  registration  of  US  and  CT  using  the  tracking 
information  (obtained  from  the  Polaris  camera).  The  registration  is  included  in  the 
data  acquisition  software,  which  enabled  fast  and  smooth  data  collection  and 
visualization.  The  application  interfaces  with  the  ultrasound  machine  and  the  tracker 
(optical  or  magnetic),  synchronizes,  records,  and  visualizes  the  data.  It  is 
implemented  in  C++  with  a  multi-threaded  scheme  for  better  performance.  We  have 
also  developed  a  simple  and  intuitive  graphical  user  interface  (GUI). 

B-mode  images  and  transducer  locations  are  continuously  stored  in  two  circular 
buffers  in  separate  threads  with  the  highest  rate  possible.  A  time-stamp  with 
microsecond  accuracy  is  also  recorded  along  with  the  captured  data.  A  constant 
delay  calculated  off-line  is  applied  to  synchronize  the  stream  of  data.  Acquisition  of 
B-mode  images  is  fast,  and  the  region  of  interest  can  be  scanned  within  a  few 
seconds.  The  maximum  rate  of  storing  data  is  bounded  with  the  minimum  of  frame 
rate  of  the  B-mode  images  and  the  tracker.  However,  the  user  can  set  the  rate  of  data 
storage  to  a  lower  rate  through  the  GUI  if  needed.  Due  to  the  high  data  rate  and  large 
frame  sizes  captured  per  second  (typically  10  to  30  frames  per  second),  immediate 
saving  to  the  hard  drive  may  result  in  unexpected  delays  or  loss  of  frames. 
Therefore,  the  data  is  first  stored  in  the  random  access  memory  (RAM)  of  the 
computer  and  stored  in  hard  disk  only  when  the  data  collection  is  stopped. 

Once  the  B-mode  scan  is  acquired,  the  program  constructs  a  3D-US  volume.  The 
volume  is  presented  as  three  orthogonal  slices  that  can  be  freely  translated  and 
rotated  in  space  as  shown  in  Figure  2.  VTK  (Kiteware  Inc.)  is  mainly  used  for  the 
visualization  and  reslicing  tasks  (Figure  3).  The  real-time  strain  images  are  captured 
in  the  same  manner  except  that  the  strain  plane  is  overlaid  on  3D-US.  In  addition, 
the  program  can  sends  commands  to  ultrasound  machine  for  capturing  sequences  of 
RF  data.  When  the  ultrasound  machine  receives  the  command,  captures  and  saves 
one  sequence  of  RF  frames,  while  our  program  records  the  tracking  information. 
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(a)  Resliced  CT  (b)  Magnified  CT  (c)  B-mode  image 


Figure  3.  Registration  of  CT  and  ultrasound,  (a)  Shows  a  CT  image,  sliced  in  3D 
(using  interpolation)  to  match  the  ultrasound  image,  whose  orientation  is  known 
from  the  tracking  information.  In  (b),  the  square  region  of  (a)  is  magnified,  (c)  shows 
the  corresponding  ultrasound  image. 

In  addition,  we  performed  some  efforts  to  improve  the  tracking-based  registration 
results  of  Figure  3  by  generating  a  pseudo-ultrasound  image  from  the  CT  volume 
and  registering  the  CT  and  ultrasound  through  the  pseudo-ultrasound,  but  the 
preliminary  results  were  not  satisfactory,  and  we  decided  to  use  the  tracking-based 
registration  results  of  Figure  3. 


d.  Compare  different  combination  of  US,  USE  and  CT  for 

delineation  (months  3-18). 

We  have  done  preliminary  studies  on  comparing  the  performance  of  the  operator 
with  different  imaging  modalities.  While  it  is  intuitive  that  adding  ultrasound  and 
elastography  enhance  the  delineation  of  the  lumpectomy  bed,  we  have  not  achieved 
a  conclusive  result  yet. 


e.  Optimize  USE  code  for  using  the  human  data  (months  3-24).  I  have 

developed  a  novel  USE  technology  and  software  under  a  Breast  Cancer  Research 
Foundation  research.  The  work  thus  far  has  been  extremely  promising:  a  manuscript 
has  already  been  accepted  for  publication  in  the  IEEE  Trans.  Med.  Imag.  I  will  be 
further  improving  my  USE  implementation  to  enhance  visualization  of  the 
lumpectomy  bed  and  ductal  tissue. 
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Extremely  promising  results  have  been  obtained  in  this  area.  We  have  developed 
new  elastography  techniques  that  generate  superb  images  from  the  human  data.  In  J1 
(cited  below  and  attached  to  the  report),  we  introduced  a  2D  strain  imaging 
technique  based  on  minimizing  a  cost  function  using  dynamic  programming  (DP). 
The  cost  function  incorporates  similarity  of  echo  amplitudes  and  displacement 
continuity.  Since  tissue  deformations  are  smooth,  the  incorporation  of  the 
smoothness  into  the  cost  function  results  in  reduced  decorrelation  noise.  As  a  result, 
the  method  generates  high  quality  strain  images  of  freehand  palpation  elastography 
with  up  to  10%  compression,  showing  that  the  method  is  more  robust  to  signal 
decorrelation  (caused  by  scatterer  motion  in  high  axial  compression  and  non-axial 
motions  of  the  probe)  in  comparison  to  the  standard  correlation  techniques.  The 
method  operates  in  less  than  1  second  and  is  thus  also  potentially  suitable  for  real 
time  elastography.  In  J2  (cited  below  and  attached  to  the  report),  we  further 
proposed  an  analytic  minimization  technique  for  strain  estimation  that  generates 
smooth  strain  images  from  the  DP  USE  images.  In  Cl  (cited  below  and  attached  to 
the  report),  we  proposed  to  exploit  tracking  information  to  select  the  best  frames  for 
USE,  and  showed  that  optimal  frame  selection  enhances  the  USE  quality.  In  C2 
(cited  below  and  attached  to  the  report),  we  introduce  a  technique  for  calculating  a 
displacement  field  from  three  frames  of  ultrasound  RF  data.  To  this  end,  we  first 
introduced  constraints  on  variations  of  the  displacement  field  with  time  using 
mechanics  of  materials.  These  constraints  are  then  used  to  generate  a  regularized 
cost  function  that  incorporates  amplitude  similarity  of  three  ultrasound  images  and 
displacement  continuity.  We  optimize  the  cost  function  in  an  expectation 
maximization  (EM)  framework.  Iteratively  reweighted  least  squares  (IRLS)  is  used 
to  minimize  the  effect  of  outliers.  We  show  that,  compared  to  using  two  images,  the 
new  algorithm  reduces  the  noise  of  the  displacement  estimation.  The  displacement 
field  is  used  to  generate  strain  images  for  quasi-static  elastography. 


Cl.  Rivaz,  H.,  Foroughi,  P.,  Fleming,  I.,  Zellars,  R.,  Boctor,  E.,  Hager,  G.,  Tracked 
Regularized  Ultrasound  Elastography  for  Targeting  Breast  Radiotherapy,  Medical 
Image  Computing  and  Computer  Assisted  Intervention,  MICCAI,  London,  UK, 
Sept.  2009,  pp  507-515  [acceptance  rate:  33%]  Conference  listed  in  PubMed  and 
papers  treated  and  cited  similar  to  high-impact  factor  journal  papers. 
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Multiple  Images”,  IEEE  Trans.  Medical  Imaging  [acceptance  rate:  30%].  Conference 
listed  in  PubMed  and  papers  treated  and  cited  similar  to  high-impact  factor  journal 
papers. 

Jl.  Rivaz,  H.,  Boctor,  E.,  Foroughi,  P.,  Zellars,  R.,  Fichtinger,  G.,  Hager,  G., 
Ultrasound  Elastography:  a  Dynamic  Programming  Approach,  IEEE  Trans.  Medical 
Imaging,  Oct.  2008,  vol.  27  pp  1373-1377 

J2.  Rivaz,  H.,  Boctor,  E.,  Choti,  M.,  Hager,  G.,  Real-Time  Regularized  Ultrasound 
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f. 

18-36). 


Performance  optimization  and  refinements  of  subsystems  (months 


The  ultrasound  elastography  subsystem  is  almost  finalized.  We  have  made  the  code 
also  available  to  the  public1,  so  that  groups  who  work  in  different  applications  of 
ultrasound  elastography  can  exploit  our  novel  and  high-quality  elastography 
estimation  technique.  We  could  confidently  say  that  the  method  can  be  implemented 
commercially  with  small  modifications. 


2.  Obtain  tracked  US  and  USE  scans  of  the  patients  weekly  during  the 
radiotherapy  (months  2-14).  This  data  will  be  acquired  from  the  same  patients 
from  whom  US  and  USE  data  were  acquired  in  step  l.b. 

As  mentioned  in  the  lb  Section,  elastography  and  ultrasound  images  are  both 
acquired  from  patients. 


3.  Develop  US  to  US  registration  (month  25-27).  This  task  is  partly  solved  in 
task  l.c.  (US  to  CT  registration).  As  mentioned,  US  to  US  registration  is  one  of  the 
steps  required  in  US  to  CT  registration.  I  will  optimize  this  step  to  work  with  breast 
images,  as  opposed  to  simulation  images  in  task  l.c. 

As  mentioned  before,  registration  of  the  pseudo-ultrasound  to  the  ultrasound  images 
are  required  to  perform  CT  to  ultrasound  registration.  Please  refer  to  the  task  l.c  for 
more  details. 


KEY  RESEARCH  ACCOMPLISHMENTS: 

•  Development  of  a  novel,  real-time,  robust  and  high  quality  elastography 
techniques 

•  Development  of  a  tracked  ultrasound  system  using  magnetic  trackers 

•  Development  of  a  technique  for  3D  volumetric  ultrasound  imaging  using  US 
image  data 

•  Study  of  ultrasound  images  of  lumpectomy  cavity 

•  Study  of  elastography  images  of  lumpectomy  cavity 

•  Data  acquisition  from  more  than  15  patients 

•  Developed  the  2D  AM  method  which  calculates  high  quality  2D  strain  images 
in  real-time 


1  Available  online  at  www.cs.ihu.edu/~rivaz/Ultrasound  Elastography 
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•  Introduced  ElastMI  (Elastography  using  Multiple  Images),  a  novel  method 
which  generates  high  quality  strain  images  by  utilizing  multiple  ultrasound  images. 


REPORTABLE  OUTCOMES: 
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Training: 


The  training  obtained  during  the  year  of  2010  was  mainly  in  the  following  forms: 

•  Reading  papers  related  to  partial  breast  radiation  therapy  techniques 

•  Helping  in  preparation  of  Institutional  Review  Boards  (IRB)  for  patient  trials 

•  Meetings  with  radiation  oncologists  on  the  problems  of  current  breast 
radiotherapy  workflow 

•  Analyzing  and  enhancing  the  CT,  ultrasound  and  strain  images  obtained  from 
the  lumpectomy 

•  Setting  up  the  image  acquisition  setup  for  imaging  lumpectomy  patients  and 
helping  in  CT  and  ultrasound  data  collection. 


CONCLUSION: 

Data  acquisition  from  lumpectomy  patients  went  smoothly.  We  acquired  the 
ultrasound  data  few  minutes  before  the  CT  scan  was  acquired.  This  means  that  we 
did  not  alter  the  current  medical  work-flow.  Ultrasound  data  acquisition  took  only 
few  minutes,  minimizing  the  cost  and  patient  discomfort.  Magnetic  tracking 
provided  ease  of  use  and  enough  reliability  for  registration  of  ultrasound  to  CT 
images  and  for  reconstructing  3D  volumes  from  2D  data.  The  novel  ultrasound 
elastography  technique  has  many  advantages  compared  to  the  previous  methods.  We 
chose  the  novel  application  of  the  lumpectomy  cavity  localization  as  the  hard  scar 
tissue  around  the  lumpectomy  is  relatively  thin  and  demands  a  high  resolution 
elastography  method.  Also,  incoherent  fluid  motions  in  the  cavity  causes  large 
decorrelations,  requiring  a  robust  method.  Due  to  the  hard  and  noisy  ultrasound 
images  acquired  from  lumpectomy  patients,  we  developed  robust  image  analysis 
techniques,  namely  the  Dynamic  Programming  (DP)  elastography  method  as 
explained  before,  the  Analytic  Minimization  (AM),  Kalman  filtering  and  Iteratively 
Reweighted  Least  Squares  (IRLS)  USE  methods  as  explained  before,  the 
Elastography  Using  Multiple  Images  (ElastMI)  framework  for  exploiting  multiple 
images  to  enhance  USE  quality,  and  the  Expectation  Maximization  (EM)  method  for 
3D  volumetric  ultrasound  generation  using  the  ultrasound  scatterer  information.  We 
also  believe  that  USE  is  useful  also  for  imaging  breast  tumors  to  reduce  the  number 
of  biopsy  referrals,  as  USE  is  an  extremely  more  convenient  and  less  expensive 
option. 
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Ultrasound  Elastography:  A  Dynamic 
Programming  Approach 

Hassan  Rivaz*,  Emad  Boctor,  Pezhman  Foroughi,  Richard  Zellars,  Gabor  Fichtinger,  and  Gregory  Hager 


Abstract — This  paper  introduces  a  2-D  strain  imaging  technique 
based  on  minimizing  a  cost  function  using  dynamic  programming 
(DP).  The  cost  function  incorporates  similarity  of  echo  ampli¬ 
tudes  and  displacement  continuity.  Since  tissue  deformations  are 
smooth,  the  incorporation  of  the  smoothness  into  the  cost  function 
results  in  reduced  decorrelation  noise.  As  a  result,  the  method 
generates  high-quality  strain  images  of  freehand  palpation  elas¬ 
tography  with  up  to  10%  compression,  showing  that  the  method 
is  more  robust  to  signal  decorrelation  (caused  by  scatterer  motion 
in  high  axial  compression  and  nonaxial  motions  of  the  probe)  in 
comparison  to  the  standard  correlation  techniques.  The  method 
operates  in  less  than  1  s  and  is  thus  also  potentially  suitable  for 
real  time  elastography. 

Index  Terms — Dynamic  programming,  freehand  ultrasound 
(US),  real  time  strain  imaging,  regularization,  ultrasound  elastog¬ 
raphy. 


I.  Introduction 

ELASTOGRAPHY,  the  computation  of  the  spatial  varia¬ 
tion  of  the  elastic  modulus  of  tissue,  is  an  emerging  med¬ 
ical  imaging  method  with  medical  applications  such  as  tumor 
detection  [1].  This  paper  focuses  on  static  elastography,  a  well- 
known  technique  that  applies  quasi-static  compression  of  tissue 
and  simultaneously  images  it  with  ultrasound.  Through  analysis 
of  the  ultrasound  images,  a  tissue  displacement  map  can  be  ob¬ 
tained  [2],  [3].  A  least  squares  technique  is  then  typically  used 
to  generate  a  low  noise  strain  estimate  from  the  displacement 
map  [2]. 

Despite  having  numerous  potential  clinical  applications,  sev¬ 
eral  practical  challenges  have  hindered  wide  application  of  static 
elastography.  First,  signal  decorrelation  between  the  precom¬ 
pression  and  postcompression  images  induces  significant  noise 
in  the  obtained  displacement  map  and  is  one  of  the  major  lim¬ 
iting  factors  in  elastography  [4].  Major  sources  of  signal  decor- 
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relation  are  scatterer  motion  in  high  axial  compression,  nonaxial 
motions  of  the  probe,  and  physiologic  motion.  Most  elastog¬ 
raphy  techniques  estimate  local  displacements  of  tissue  based 
on  correlation  analysis  of  radio-frequency  (RF)  echoes  [2],  [3]. 
Large  windows  are  required  to  reduce  the  variance  (i.e.,  noise) 
of  the  estimated  displacement  and  to  avoid  ambiguity  in  time 
delay  estimation,  especially  when  tracking  a  motion  that  ex¬ 
ceeds  one  wavelength.  At  the  same  time,  signal  decorrelation 
within  large  windows  limits  the  tolerable  level  of  compression 
[2]— [4] .  To  reduce  signal  decorrelation,  stretching  methods  have 
been  proposed  [5],  [6],  which  are  computationally  expensive 
and  are  not  suitable  for  real-time  elastography.  Moreover,  large 
errors  due  to  false  peaks  and  smaller  errors  due  to  jitter  [7]  limit 
the  performance  of  correlation  techniques. 

Second,  in  many  methods,  the  compression  is  applied  by  a 
mechanical  actuator  in  order  to  generate  an  excitation  that  mini¬ 
mizes  signal  decorrelation  [1],  [8]  or  because  accurate  motion  is 
otherwise  required  by  the  particular  elastography  technique  [9]. 
Freehand  palpation  elastography  is  a  much  more  attractive  alter¬ 
native,  as  it  requires  no  extra  hardware  and  provides  ease  of  use. 
It  has  attracted  increasing  interest  in  recent  years  [8],  [10]— [13], 
however  it  introduces  additional  sources  of  signal  decorrelation 
caused  by  operator’ s  hand  unwanted  motion. 

Third,  elastography  is  computationally  expensive,  making  it 
challenging  to  display  elastograms  in  real  time.  Real-time  elas¬ 
tography  provides  the  feedback  to  the  operator  to  best  capture 
the  region  of  interest  in  the  elastogram  and  is  required  for  image 
guided  surgical  operations  that  can  potentially  use  elastograms. 
Combined  autocorrelation  method  [14]  and  phase  zero  estima¬ 
tion  [15]  are  the  first  work  that  generate  real-time  elastograms. 
Hall  et  al.  [12]  have  presented  a  real-time  elastography  system 
where  tissue  compression  is  performed  by  freehand  palpation 
based  on  a  2-D  block  matching  algorithm.  Dynamic  program¬ 
ming  is  used  for  one  A-line  of  the  image  for  guiding  the  block 
matching  algorithm  [16].  While  these  methods  use  the  displace¬ 
ment  of  each  window  to  confine  the  search  range  for  the  neigh¬ 
boring  windows,  the  displacement  of  each  window  is  calculated 
independently  and  hence  are  sensitive  to  signal  decorrelation. 

In  work  closely  related  to  this  paper,  Pellot-Barakat  et  al.  [17] 
have  proposed  minimizing  an  energy  function  that  combines 
constraints  of  conservation  of  echo  amplitude  and  displacement 
continuity.  Since  data  alone  can  be  insufficient  to  solve  ambi¬ 
guities  due  to  signal  decorrelation,  the  physical  priors  of  tissue 
motion  continuity  increases  the  robustness  of  the  technique.  The 
RF  data  is  first  upsampled  by  a  factor  of  four  in  the  axial  direc¬ 
tion.  The  image  is  then  subdivided  into  four  parts  and  a  coarse 
displacement  map  is  calculated  for  each  part  iteratively.  Each 
part  is  subsequently  divided  into  four  parts  and  the  displace- 
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merit  of  each  part  is  calculated  by  the  same  iterative  technique 
using  the  displacement  of  the  parent  grid  as  an  initial  guess.  The 
method  is  shown  to  generate  accurate  low  noise  displacement 
fields.  However,  the  computation  time  is  reported  to  be  more 
than  1  min  for  a  strain  image  that  is  less  than  half  of  the  number 
of  pixels  in  the  strain  images  generated  in  this  paper.  Hence,  the 
method  is  not  immediately  suitable  for  real  time  elastography. 

The  contribution  of  this  paper  is  the  demonstration  of  the 
feasibility  of  an  elastography  technique  based  on  dynamic  pro¬ 
gramming  (DP)  for  image  matching  [18].  Compared  to  other 
optimization  techniques,  DP  is  an  efficient  noniterative  method 
of  global  optimization  [19],  [20].  However,  it  can  only  be  used 
to  optimize  causal  cost  functions  (Section  II). 


d  (displacement) 
-10  12  3  4 


Fig.  1.  I  n  the  I  eft,  values  of  g(i)  and  gr(t+d)  corresponding  to  precompression 
and  postcompression  RF  data  are  compared.  Right  shows  the  cost  function  C 
of  (4)  (white  and  black  represent  low  and  high  cost  values,  respectively). 


II.  One-Dimensional  Displacement  Estimation  Using  DP 

Devising  a  DP  algorithm  for  optimization  involves  the  fol¬ 
lowing. 

1)  Breaking  the  total  optimization  cost  into  a  sum  of  indi¬ 
vidual  costs,  such  that  each  cost  corresponds  to  a  discrete 
decision.  The  decisions  should  follow  each  other  sequen¬ 
tially  and  the  cost  corresponding  to  each  decision  should 
only  depend  on  the  previous  and  not  the  future  decisions 
(causality). 

2)  Determining  what  decisions  are  possible  at  each  stage. 

3)  Writing  a  recursion  on  the  optimal  cost  from  the  first  stage 
to  the  final  stage. 

We  first  consider  the  problem  of  1-D  strain  estimation  with 
1-D  smoothness  regularization.  Consider  two  echo  signals  g(i) 
and  g'(i)  correspondi  ng  to  two  A  -I  i  nes  acqui  red  before  and  after 
compression  (Fig.  1,  left),  each  signal  sampled  at*  =  1, 

The  difference  between  the  two  signals  A  can  be  quantified 
using  sum  of  absolute  differences  (SAD),  which  is  computation¬ 
ally  inexpensive  to  compute  and  has  been  shown  to  have  good 
robustness  against  outliers  [21],  [20] 

A(i,d)  =  \g(i)  -g'(i  +  d)\  (1) 

where  dmin  <  d  <  dmax  is  the  displacement  at  the  sample  * 
(Fig.  1,  left)  and  dmm  and  dmaK  specify  the  allowed  displace¬ 
ment.  Gains  of  RF  data  can  be  changed  in  ultrasound  machines 
to  improve  visualization.  To  reduce  the  effect  of  these  changes 
on  A,  both  precompression  and  postcompression  ultrasound  im¬ 
ages  are  divided  by  the  maximum  value  of  one  of  the  images. 
The  smoothness  of  the  displacements  is  S 

S(di,di-1)  =  (di-di-1)k  (2) 

where  di  is  the  displacement  at  the  sample  *  and  di-i  is  the 
displacement  at  the  sample  *  -  1  of  the  g(i).  To  avoid  large 
jumps  in  the  displacement,  S  should  be  strictly  convex 

q  ( dn  -  di-i)k  +  (1  -  a)  ( di2  -  di-i)k  > 

[adn  +  (1  -  a)di2  -  di-i]k  ,  0<a<l  (3) 

i.e.,  a  small  jump  and  a  large  jump  (left-hand  side)  are  penalized 
more  than  two  medium  jumps  (right-hand  side).  This  holds  for 
even  k,  we  choose  k  —  2:  for  k  >  2  larger  jumps  are  more 
heavi  I  y  penal  ized  which  adversel  y  affects  contrast  to  noi  se  rati  o. 


The  cost  function  Cat  a  point*  and  associated  displacement^ 
is  defined  as  a  recursive  function 

C(*,  di)  =  min  {C(*  -  1,  di- 1)  +  wS(di,  di- 1)}  +  A (*,  di) 

di- 1 

(4) 

where  w  is  a  regularization  weight  which  governs  smooth¬ 
ness.  The  study  of  its  effect  on  the  estimated  displacement  is 
postponed  to  the  discussion  of  2-D  displacement  estimation 
in  Section  1 1 1 -A .  The  values  of  the  C  function  are  stored  in  a 
(rfmax  m  i  ]  i  +  1)  x  m  matrix  (Fig.  1,  right). 

G  eneral  ly,  the  opti  mum  val  ue  of  di- 1  shoul d  be  sought  i  n  the 
entire  [dmin,  dmax]  range.  However,  si  nee  the  strain  value  is  low 
in  elastography,  it  is  expected  and  desi  red  that  at  each  sample  of 
RF  data,  the  change  between  the  displacement  of  a  sample  and 
its  previous  sample  is  not  more  than  1.  Therefore,  the  search 
range  of  opti  mum  val  ue  for  dt-i  i  s  I  i  mi  ted  to  the  three  val  ues  of 
di- 1,  di  and  di  +  1,  which  results  in  a  significant  gain  in  speed. 
This  limit  on  the  search  range  does  not  affect  the  results  even  in 
a  high  strain  of  10%:  Ad  is  zero  for  nine  samples  and  one  for 
the  tenth  sample  on  average.  The  value  of  di- 1  that  minimizes 
(4)  is  also  "memoized"  [19]  in  a  function  M  for  later  use 

M(i,  di)  =  arg min  {C(*  -  1,  di- 1)  +  wS(di,  .  (5) 

di- 1 

The  cost  function  C  is  calculated  for*  =  1  •  •  -m.  The  minimum 
cost  at  *  =  m  gives  the  displacement  of  this  point,  which  is 
traced  back  to  *  =  1  using  the  M  function  to  calculate  all  the 
displacements  (D) 

D(i)  =argmin{C'(*,(fi)} ,  i  —  m 

di 

D(i)  —M(i  +  1,£>(*  +  1)),  i  —  —  1.  (6) 

Thedisplacementmapof  all  A -lines  is  calculated  using  the  same 
procedure  independently.  I  n  Section  1 1 1 ,  we  present  a  method  for 
coupling  adjacent  A-lines. 

A.  Hierarchical  Search  and  Subpixel  D isplacement E stimation 

Further  speedup  is  achieved  by  downsampling  the  signal  <?(*) 
by  a  factor  of  (3  to  <?*(*),  and  comparing  it  with  the  unaltered 
signal  g'{i).  This  is  done  by  simply  skipping  (3  -  1  samples 
from  g(i)  and  performing  DP  on  the /3th  sample,  as  illustrated 
in  Fig.  2  left.  This  gen erates  integer  displacement esti m ati o n s  at 
m//3  samples.  The  displacement  of  the  skipped  samples  is  then 
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Fig.  2.  I  n  the  left,  the  cost  function  C  is  shown  when  DP  is  performed  on  g*(i) 
[g(i]  downsampled  by  a  factor  of  0)  and  gr(i)  (not  downsampled).  Hashed 
squares  indicate  no  cost  calculation  is  performed  due  to  downsampling  of  g(i), 
and  white  and  black  representing  low  and  high  cost  values  respectively.  Dis¬ 
placement  is  calculated  at  m/ (3  samples  in  this  stage  (p  =  3  in  this  figure). 
In  right,  a  new  cost  function  around  the  optimum  path  of  the  first  stage's  cost 
function  (the  dashed  line)  is  created,  giving  a  l/j  =  1/2  pixel  displacement 
accuracy  at  m  samples. 


simply  approximated  by  the  linear  interpolation  of  two  neigh¬ 
boring  points  whose  displacements  are  calculated,  as  an  initial 
guess  for  the  next  step. 

The  displacement  estimates  are  then  refined  to  subpixel  dis¬ 
placement  estimation  at  all  m  samples.  The  original  signal  g(i ) 
(not  downsampled)  is  compared  with  g'{i  +  d)  upsampled  by  a 
facto r  of -y  (Fig.  2  right)  using  parabolic  interpolation.  R epeati n g 
the  refinement  procedure  n  times  results  in  a  refinement  factor 
of  1/7". 

I  n  cross  correlation  methods,  subsample  displacement  is  usu¬ 
ally  achieved  by  interpolation  of  the  correlation  function  [22], 
which  is  subject  to  bias  and  jitter  [22],  [23],  Here,  we  interpo¬ 
late  the  original  RF  data  instead,  which  is  shown  to  have  similar 
performance  [23],  Although  cosine-fit  outperforms  parabolic-fit 
interpolation  in  terms  of  bias  and  jitter  [22],  [23],  the  latter  is 
used  here  for  computational  simplicity. 


strain  %  width  (mm) 

<e)  (f) 

Fig.  3.  (a)-(c)  strain  images  obtained  from  freehand  palpation  of  the  phantom 
using  cross  correlation,  cross  correlation  with  a  3  x  3  median  filter  applied  on 
the  displacement  image  and  1-D  DP  respectively.  The  target  window  is  fixed  on 
the  lesion  and  the  background  window  is  moved  to  allow  multiple  CNR  calcu¬ 
lation.  (d)  Normalized  CNR  values  of  the  lesion,  obtained  by  dividing  each  bin 
by  the  total  of  36  C  N  R  measurements,  (e)  SN  R  values  of  the  cross  correlation 
and  1-D  DP  techniques,  (f)  Strain  images  obtained  from  freehand  palpation  of 
the  phantom  using  2-D  DP. 


B.  Results 

For  experimental  evaluation,  RF  data  was  acquired  from 
an  Antares  Siemens  system  (Issaquah,  WA)  with  a  7.27-M  Hz 
linear  array  at  a  sampling  rate  of  40  M  Hz.  For  the  purposes  of 
comparison,  strain  images  were  also  calculated  using  a  standard 
cross  correlation  method  with  a  3-mm  window  size  and  80% 
overlap  and  a  three  point  parabolic  interpolation  to  find  the  sub¬ 
sample  location  of  the  correlation  peak  [22],  Linear  regression 
with  a  5-sample  window  is  performed  on  the  displacement  field 
to  calculate  strain.  Normalization  was  performed  to  decrease 
the  dynamic  range  of  the  strain  images:  any  strain  value  outside 
s  =L  3(7  was  set  to  s  ±  3a  to  eliminate  the  outliers  in  the  strain 
map  (s  and  a  are  the  mean  and  standard  deviation  of  the  strain 
values  across  the  whole  image).  The  unitless  performance 
metric  signal-to-noise  ratio  (SNR)  and  contrast-to-noise  ratio 
(CNR)  were  calculated  according  to  [2] 


(7) 


where  st  and  sb  are  the  spatial  strain  average  of  the  target  and 
background,  of  and  of  are  the  spatial  strain  variance  of  the 


target  and  background,  and  s  and  a  are  the  spatial  average  and 
variance  of  a  window  in  the  strain  image,  respectively. 

In  the  first  experiment,  a  breast  el  astography  phantom  (CIRS, 
Norfolk,  VA)  with  a  lesion  of  10  mm  diameter  and  three  times 
stifferthan  the  background  was  palpated  freehand.  In  consecu¬ 
tive  images,  where  axial  compression  is  low  and  there  is  little 
nonaxial  motion,  both  methods  perform  well.  However,  as  the 
axial  compression  and  attendant  nonaxial  motion  increase,  the 
DP  method  outperforms  the  cross  correlation  method.  Fig.  3(a) 
shows  the  strain  image  obtained  with  cross  correlation.  In 
Fig.  3(b),  a  3  x  3  median  filter  is  applied  to  the  displacement 
measurements,  before  differentiation,  as  a  2-D  continuity 
check.  Fig.  3(c)  shows  the  strain  image  obtained  with  the  1-D 
DP  method.  A  high  level  of  lateral  motion,  slightly  more  than  2 
A-lines,  at  the  left  of  the  image  and  high  axial  strain  cause  the 
cross  correlation  method  to  fail.  To  calculate  the  CNR  values 
the  target  window  was  selected  as  specified  in  the  figure.  The 
background  window  was  then  moved  across  the  strain  image 
(with  3.8  mm  margin  from  all  four  sides  and  from  the  lesion 
where  the  strain  is  expected  to  vary  considerably)  to  allow 
for  a  more  comprehensive  CNR  measurement.  The  histogram 
of  Fig.  3(d)  shows  that  1-D  DP  gives  better  CNR  values:  the 
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mean  value  of  the  36  CNR  measurements  for  cross  correlation, 
cross  correlation  with  the  3  x  3  median  filter,  and  1-D  DP  are, 
respectively,  2.60,  3.98,  and  6.24.  The  standard  deviation  value 
of  CNR  for  cross  correlation,  cross  correlation  with  the  3  x  3 
median  filter,  and  1-D  DP  are,  respectively,  2.08,  2.70,  and 
4.27,  reflecting  the  changes  in  strain  across  the  image  caused 
by  medium  inhomogeneity  and  nonuniform  loading  condition. 

To  obtain  a  strain  filter  [2],  a  CIRS  elasticity  QA  phantom 
with  the  Young's  modulus  of  33  kPa  was  compressed  in  24 
steps,  each  step  0.005  in.  The  experiment  was  performed  far 
from  the  lesions  of  the  phantom  to  generate  close  to  uniform 
strain  due  to  a  uniform  compression.  The  strain  map  between 
the  first  frame  and  all  other  frame  was  calculated  using  the  cross 
correlation  and  DP  methods.  The  SNR  metric  was  calculated  in 
a  small  window  located  at  the  top  center  of  the  image,  where 
strain  is  approximately  constant.  Fig.  3(e)  shows  that  the  1-D 
DP  method  has  a  higher  dynamic  range,  an  important  elastog- 
raphy  performance  metric  [2], 


Fig.  4.  Two-dimensional  DP  results  of  freehand  palpation  of  the  breast 
phantom,  (a)  Strain  image  with  the  target  window  and  four  background 
windows  for  CNR  calculation.  (b)CNR  values  between  the  target  window  and 
the  four  background  windows  on  the  top,  right,  bottom  and  left  of  the  lesion 
calculated  for  different  regularization  coefficient  values  w  in  (10). 


Cj+i(da,  di,  i)  and  is  deleted  from  the  memory  afterward.  This 
makes  the  amount  of  memory  required  to  store  the  cost  function 
values  independent  of  the  number  of  A-lines. 


III.  2-D  Displacement  Estimation 

Until  now,  we  have  assumed  pure  axial  compression  inde¬ 
pendently  estimated  on  each  A-line.  Flowever,  lateral  displace¬ 
ment  in  a  soft  material  is  inevitable  even  when  it  undergoes  pure 
axial  compression.  This  displacement  is  related  to  the  Poisson's 
ratio,  which  describes  the  material  compressibility.  Also,  free¬ 
hand  palpation  is  rarely  a  pure  compression  and  thus  also  re¬ 
sults  in  nonaxial  tissue  motion.  Asa  result,  a  2-D  smoothness 
regularization  that  considers  the  displacements  between  adja¬ 
cent  A-lines  is  more  natural.  The  DP  algorithm  of  Section  II  is 
modified  hereto  allow  for  2-D  displacement  estimation  and  2-D 
smoothness  regularization. 

Assuming  that  ultrasound  images  consist  of  n  A-lines,  the 
distance  between  the  pre  and  postcompression  signals  is 

A (i,j,da,di)  =  \gj(i)  -  g'j+dl(i  +  da)\  (8) 

where  da^mm  £  da  £  da^ max  3nd  c^min  ^  di  dk max  ^re  the 

axial  and  lateral  displacements,  respectively,  and  j  =  1  ■■■n 
refers  to  jth  A-line  and  i  =  1  •  •  -m 

S (dai ,  dh ,  daj_1 ,  dh_t )  -  (dai  - da i_l)1  +  (du  - dk_ L )2  (9) 

is  the  smoothness  regularization  with  subscripts  a  and  /referring 
to  axial  and  lateral.  The  cost  function  atthefth  sample  of  the  jth 
A-line  is 


Cj(da,di,i)  =  A(da,di,i)  +  min 


i&l 


Cj{6a,6u%-\)  +  Cj^{6a,6u%) 


wS(da,di,Sa,Si ) 


(10) 


For  memoization,  6a  and  bt  values  that  minimize  the  cost 
function  are  stored  for  all  da,  dt  and  i  values.  The  specific  form 
of  the  cost  function  allows  the  calculation  of  the  displacement 
of  each  A-line  using  the  cost  values  of  the  previous  A-line.  The 
cost  function  of  the  jth  line,  Cj(da,di,i),  is  calculated  and  is 
minimized,  resulting  in  its  displacement  map.  The  Cj(da,di,  i) 
function  is  also  used  for  the  calculation  of  the  next  cost  function 


A.  Results 

To  study  the  effect  of  the  regularization  weights  on  theCN  R, 
the  breast  phantom  is  palpated  freehand.  For  two  RF  frames, 
theelastogram  is  obtained  using  the  2-D  DP  algorithm  with  dif¬ 
ferent  w  values.  CNR  is  calculated  between  the  shown  target 
window  and  thefour  background  windows  on  top,  right,  bottom, 
and  left  of  the  lesion  (Fig.  4).  At  low  w  values,  CNR  is  low 
because  of  high  noise  in  the  windows,  while  at  high  w  values 
CNR  drops  because  high  displacement  changes  are  heavily  pe¬ 
nalized.  3  <  w  <  50  is  optimizing  the  tradeoff  between  noise 
and  contrast  to  maximize  CNR,  both  in  the  lateral  (background 
windows  on  the  right  and  left  of  the  target)  and  axial  (back¬ 
ground  windows  on  the  top  and  bottom  of  the  target)  directions. 
Since  the  background  windows  are  close  to  the  lesion,  the  strain 
within  each  window  is  not  expected  to  be  constant  even  though 
the  phantom  is  homogeneous  within  them.  This  variation  in  the 
ground  truth  strain  will  be  reflected  as  noise  in  the  CNR  cal¬ 
culation  which  is  undesirable.  Flowever,  we  have  selected  the 
windows  such  that  they  best  capture  the  effect  of  w  on  CNR 
close  to  an  inhomogeneity,  which  also  contains  some  resolution 
information. 

Fig.  3(f)  shows  the  strain  image  obtained  using  the  2-D  DP 
method  using  the  same  two  frames  that  are  used  to  generate  the 
strain  images  in  Fig.  3(a)-(c).  The  CNR  values  [Fig.  3(d)]  are 
calculated  for  the  same  target  and  36  background  windows  as 
before,  giving  a  mean  of  8.96  and  a  standard  deviation  of  5.75. 
Since  the  elasticity  QA  phantom  cannot  be  compressed  more 
than  4%,  we  use  the  breast  phantom  for  experimental  evalua¬ 
tion  of  the  strain  filter  of  the  2-D  DP  method.  Fig.  3(e)  shows 
the  SN  R  values,  showing  no  degradation  of  the  SN  R  even  at  a 
high  strain  of  10%.  Comparing  these  results  with  the  1-D  DP 
and  cross  correlation  results,  a  significant  increase  in  the  image 
quality,  CNR  value,  and  maximum  allowed  strain  is  achieved. 

Substituting  other  computationally  moreexpensivesimilarity 
measures  like  normalized  cross  correlation  in  the  A  function  re¬ 
sulted  in  no  significant  difference  in  the  performance.  Currently, 
the  algorithm  takes  0.72  s  to  calculate  the  displacement  map  of 
each  pixel  in  an  image  with  1000  x  100  pixels  with  maximum 
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axial  displacement  of  10  samples  (1%  strain)  and  maximum  lat¬ 
eral  displacement  of  ±1  A-lines  on  a  3.8  G  Hz  P4  C  PU .  The  cur¬ 
rent  implementation  is  in  M  ATLAB  with  the  DP  optimization 
in  mex  functions. 

IV.  Discussion  and  Conclusion 

Thelateral  search  is  performed  in  the  2- D  DP  method  onlyto 
decrease  the  noi  se  and  i  ncrease  the  robustness  of  the  axi  al  strai  n : 
the  lateral  displacements  are  integer  values  and  are  not  suitable 
for  calculating  lateral  strain.  Results  of  Fig.  3  show  that  DP  is 
more  robust  to  the  signal  decorrelation  (caused  by  scatterer  mo¬ 
tion  in  high  axial  compression,  and  lateral  and  out  of  plane  mo¬ 
tions  of  the  probe)  than  standard  cross  correlation  techniques. 
This  i  s  cri  ti  cal  I  y  i  mportant  si  nee  i  ttol  erates  higher  axial  compres¬ 
sion  [Fig.  3(e)],  increasing  the  dynamic  range  of  theelastogram 
which  is  crucial  for  lesion  detection.  Nonlinear  elastic  properties 
of  tissue  also  only  appear  at  high  strain  values  [3].  It  also  gener¬ 
ates  low-noise  elastograms  using  almost  any  two  frames  in  free¬ 
hand  pal  pati  on,  given  that  they  both  bel  ong  to  the  same  compres- 
sion  or  relaxation  cycle  of  the  palpation  excitation.  Finally,  no 
postprocessing  step  such  as  median  filtering  is  required. 

The  CNR  and  SNR  metrics  seem  to  indicate  that  the  regu¬ 
larization  creates  smooth  elastograms  while  preserving  contrast. 
The  only  tunable  parameter  of  the  method,  w  in  (10),  was  kept 
constantatu?  =  10  throughoutthis  work.  Itcan  also  bevaried  be¬ 
tween  5  and  50,  as  F  i  g.  4  i  ndi  cates,  w  i  th  al  most  no  effect  on  CNR. 
This  might  indicate  that  the  w  value  optimized  for  phantom  will 
work  well  for  real  tissue.  These  features  of  the  DP,  along  with 
its  high  speed  make  it  a  promising  elastography  method. 

DP  strain  images  in  Figs.  3  and  4  show  some  stress  concen¬ 
tration  around  the  lesion  which  is  not  seen  in  the  corresponding 
cross  correlation  images.  We  are  not  sure  yet  whether  this  is  an 
artifact  or  high  strains  are  created  just  around  the  lesion  because 
of  nonlinear  mechanical  properties  of  the  phantom.  We  are  plan¬ 
ning  for  validation  of  the  estimated  displacement  of  DP  using 
simulation  and  laboratory  experiments  for  clarification.  High 
strain  is  also  seen  on  the  top  edges  in  both  cross  correlation  and 
D  P  images.  The  curved  shape  of  the  breast  phantom  is  probably 
the  reason  for  this  high  strai  n:  i  n  order  for  the  edges  of  the  probe 
to  touch  the  phantom,  the  part  of  the  phantom  just  under  middle 
of  the  probe  has  to  compress  considerably.  If  the  phantom  ma¬ 
terial  hardens  under  high  strains,  the  phantom  around  the  edges 
experiences  higher  strain.  The  absence  of  this  stress  concen¬ 
tration  in  our  experiments  with  noncurved  phantoms  seems  to 
prove  this. 

We  have  chosen  to  use  the  cross  correlation  as  a  compara¬ 
tive  benchmark  to  assess  the  potential  of  DP.  This  is  because 
cross  correlation  is  the  most  commonly  used  method  and  has 
been  shown  to  accomplish  at  least  as  accurate  results  as  any 
other  method,  and  thus  it  represents  a  “gold-standard”  [24],  [23]. 
However,  we  are  planning  for  a  comprehensive  comparison  of 
the  DP  with  other  strain  imaging  techniques.  Further  work  is  re¬ 
quired  to  study  the  effect  of  regularization  on  resolution  [25].  To 
achieve  real-time  performance  in  freehand  palpation  imaging, 
an  adaptive  search  range  selection  can  be  implemented  by  using 
the  conti  nuity  of  di  spl  acement  i  n  ti  me  to  confine  the  search.  T  he 
2-D  algorithm  can  be  extended  to  2-D  +  t  to  exploit  the  cost 
function  in  previous  time,  optimize  frame  selection  [26],  and 
incorporate  a  2-D  +t  regularization. 
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Real-Time  Regularized  Ultrasound  Elastography 
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A  bstract-  This  paper  introduces  two  real-time  elastography 
techniques  based  on  analytic  minimization  (AM)  of  regularized 
cost  functions.  The  first  method  (ID  AM)  produces  axial  strain 
and  integer  lateral  displacement,  while  the  second  method  (2D 
AM)  produces  both  axial  and  lateral  strains.  The  cost  functions 
incorporate  similarity  of  radio-frequency  (RF)  data  intensity 
and  displacement  continuity,  making  both  AM  methods  robust  to 
small  decorrelations  present  throughout  the  image.  We  also  exploit 
techniques  from  robust  statistics  to  make  the  methods  resistant  to 
large  local  decorrelations.  We  further  introduce  Kalman  filtering 
for  calculating  the  strain  field  from  the  displacement  field  given 
by  the  AM  methods.  Simulation  and  phantom  experiments  show 
that  both  methods  generate  strain  images  with  high  SNR,  CNR 
and  resolution.  Both  methods  work  for  strains  as  high  as  10%  and 
run  in  real-time.  We  also  present  in  vivo  patient  trials  of  ablation 
monitoring.  An  implementation  of  the  2D  AM  method  as  well  as 
phantom  and  clinical  RF-data  can  be  downloaded. 

Index  Terms—  Kalman  filter,  radio-frequency  (RF)  ablation, 
real-time  ultrasound  elastography,  regulariation,  robust  estima¬ 
tion,  two-dimensional  (2D)  strain. 


I.  Introduction 

ELASTOGRAPHY  involves  imaging  the  mechanical 
properties  of  tissue  and  has  numerous  clinical  appli¬ 
cations.  Among  many  variations  of  ultrasound  elastography 
[l]-[4],  our  work  focuses  on  real-time  static  elastography,  a 
well-known  technique  that  applies  quasi-static  compression  of 
tissue  and  simultaneously  images  it  with  ultrasound.  Within 
many  techniques  proposed  for  static  elastography,  we  focus  on 
freehand  palpation  elasticity  imaging  which  involves  deforming 
the  tissue  by  simply  pressing  the  ultrasound  probe  against  it. 
It  requires  no  extra  hardware,  provides  ease  of  use  and  has 
attracted  increasing  interest  in  recent  years  [5]— [10] .  Real-time 
elastography  is  of  key  importance  in  many  diagnosis  applica¬ 
tions  [11],  [6],  [12],  [8],  [13]  and  in  guidance/monitoring  of 
surgical  operations  [14]— [16]. 

Global  and  local  decorrelation  between  the  pre-  and  post¬ 
compression  ultrasound  images  compromises  the  quality  of  the 
elasticity  images.  The  main  sources  of  global  decorrelation  in 
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freehand  palpation  elastography  are  change  of  speckle  appear¬ 
ance  due  to  scatterer  motion  and  out-of-plane  motion  of  the 
probe  (axial,  lateral  and  out-of-plane  directions  are  specified  in 
Fig.  1).  Examples  of  local  decorrelation  are:  1)  a  decrease  in 
the  ultrasonic  signal  to  noise  ratio  with  depth,  2)  low  correla¬ 
tion  close  to  arteries  due  to  complex  motion  and  inside  blood 
vessels  due  to  blood  motion,  3)  extremely  low  correlation  in  le¬ 
sions  that  contain  liquid  due  to  the  incoherent  fluid  motion  [17], 
[8],  and  4)  out-of-plane  motion  of  movable  structures  within  the 
image  [17]. 

Most  elastography  techniques  estimate  local  displacements 
of  tissue  based  on  amplitude  correlation  [18],  [2]  or  phase  corre¬ 
lation  of  the  radio-frequency  (RF)  echoes  [19]— [21].  Assuming 
a  stationary  signal  model  for  the  RF  data,  the  use  of  large  cor¬ 
relation  windows  helps  to  reduce  jitter  errors  (variance)  for  all 
motion  field  estimation  techniques  studied  in  [18]  and  [22].  This 
is  intuitive  as  larger  windows  contain  more  information.  How¬ 
ever,  in  practice  RF  data  is  not  stationary  and,  for  large  deforma¬ 
tions,  the  decorrelation  increases  with  window  size.  Therefore, 
in  addition  to  reducing  the  spatial  resolution  [23],  larger  win¬ 
dows  result  in  significant  signal  decorrelation  [24],  [23],  [18]. 
Coarse-to-fine  hierarchical  search  is  used  in  [23]  to  combine  the 
accuracy  of  large  windows  with  the  good  spatial  resolution  of 
small  window.  However,  the  issue  of  signal  decorrelation  within 
the  window  remains  unresolved  in  this  approach  and  can  cause 
the  highest  level  of  the  hierarchical  search  to  fail. 

All  of  the  aforementioned  methods  either  do  not  calculate 
the  lateral  displacement  or  they  just  calculate  an  approximate 
integer  lateral  displacement.  A  2D  displacement  field  is  re¬ 
quired  to  calculate  the  thermal  expansion,  lateral  and  shear 
strain  fields  [25]  (i.e.,  reconstruct  the  strain  tensor),  Poisson’s 
ratio  and  Young’s  modulus  [26],  [27].  The  axial  resolutions  of 
ultrasound  is  determined  by  the  pulse  length,  and  the  lateral 
resolutions  is  dictated  by  the  center  frequency  of  the  excitation 
and  the  transducer  pitch.  Therefore,  the  lateral  resolution  is  of 
order  of  magnitude  lower  than  axial  resolution.  As  a  result, 
few  2D  elastography  techniques  have  been  proposed  to  date. 
Initially,  2D  motion  estimation  started  in  the  field  of  blood  flow 
estimation  using  speckle  tracking  [28].  Designed  for  blood  flow 
estimation,  these  techniques  are  not  immediately  suitable  for 
elastography  which  involves  tissue  deformation. 

Attaching  a  coordinate  system  to  the  ultrasound  probe  as  in 
Fig.  1  ,z,x,  and  y  in  the  ultrasound  image  are  generally  defined 
as  axial,  lateral  and  out-of-plane  directions.  Assume  that  the  ap¬ 
plied  compression  to  the  tissue  is  the  Z  direction,  and  attach 
a  coordinate  system  X,  Y,  Z  to  the  applied  force.  Letting  dz 
and  dw  be  the  displacements  in  the  Z  and  N  directions  where 
X_LZ,  axial  and  transverse  strains  are  ddz/dZ  and  ddjy/ON. 
In  most  experimental  setups  (including  freehand  palpation  elas¬ 
tography),  2  and  Z  are  parallel  and  N  will  be  either  lateral  or 
out-of-plane,  and  therefore  d n  cannot  be  estimated  accurately. 
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Fig.  1.  Axial,  lateral,  and  out-of-plane  directions.  The  coordinate  system  is  attached  to  the  ultrasound  probe.  The  sample  (i,j)  marked  by  x  moved  by  hty)- 
and  lisi  are,  respectively,  axial  and  lateral  displacements  and  initially  are  integer  in  DP. 


To  calculate  an  accurate  transverse  strain,  Z  and  ^  are  perpen¬ 
dicular  in  [29]  by  applying  the  compression  force  perpendicular 
to  the  ultrasound  imaging  axis.  Therefore,  transverse  strain  is  in 
the  z  direction  of  the  ultrasound  probe  and  hence  can  be  mea¬ 
sured  with  high  accuracy.  However,  such  an  experimental  setup 
is  not  possible  in  many  medical  applications.  Beam  steering  has 
been  used  to  solve  this  issue  [30].  In  freehand  palpation  elas- 
tography,  beam  steering  causes  z  and  Z  to  be  unparallel,  so 
that  a  component  of  the  dx  is  in  the  2  direction.  The  steering 
angle  determines  the  angle  between  £  and  Z.  Unfortunately, 
large  steering  angles  are  required  to  obtain  accurate  estimates 
of  lateral  strain,  which  is  possible  in  phased  arrays  and  not  in 
linear  arrays. 

Lateral  strain1  estimation  is  obtained  in  [31]  by  iteratively  cal¬ 
culating  axial  strain,  companding  RF  data  and  interpolating  in 
the  lateral  direction.  In  another  work  [32],  tissue  deformation 
is  modeled  by  locally  affine  transformations  to  obtain  both  lat¬ 
eral  and  axial  strains.  Change  of  speckle  appearance  is  taken 
into  account  by  proposing  a  Lagrangian  speckle  model  [33].  Al¬ 
though  they  provide  high  quality  lateral  strain,  these  techniques 
are  computationally  expensive  and  are  not  suitable  for  real-time 
implementation. 

Axial  strain  is  used  in  [34]  to  enhance  the  quality  of  lateral 
displacement  estimation.  Tissue  is  assumed  to  be  incompress¬ 
ible  and  isotropic  and  therefore  axial,  lateral  and  out-of-plane 
strains  should  add  to  zero.  However,  many  tissues  cannot  be 
considered  incompressible.  In  fact,  some  research  has  even  fo¬ 
cused  on  imaging  the  ratio  of  the  axial  and  lateral  strain  (i.e.,  the 
Poisson’s  ratio  iZ)  [31]. 

While  most  previously  mentioned  methods  use  tissue  mo¬ 
tion  continuity  to  confine  the  search  range  for  the  neighboring 
windows,  the  displacement  of  each  window  is  calculated  inde¬ 
pendently  and  hence  is  sensitive  to  signal  decorrelation.  Since 
data  alone  can  be  insufficient  due  to  signal  decorrelation,  Pellot- 
Barakat  et  al.  [35]  proposed  minimizing  a  regularized  energy 
function  that  combines  constraints  of  conservation  of  echo  am¬ 
plitude  and  displacement  continuity.  In  another  work  [36],  both 
signal  shift  and  scale  are  found  through  minimization  of  a  regu- 

xWe  hereafter  assume  the  applied  force  is  in  the  2  direction  (i.e.,  Z  and  z  are 
parallel)  and  therefore  we  use  the  term  lateral  strain  instead  of  the  term  trans¬ 
verse  strain. 


larized  cost  function.  The  computation  time  of  these  methods 
is  reported  to  be  few  minutes  and  therefore  are  not  immedi¬ 
ately  suitable  for  real  time  elastography.  In  [37]  and  [38],  few 
phase-based  methods  are  regularized  and  strain  and  elasticity 
modulus  images  are  obtained.  The  regularization  term  is  the 
Laplacian  (second  derivative)  of  the  motion  field  and  is  spa¬ 
tially  variant  based  on  the  peak- value  of  the  correlation  coef¬ 
ficient.  The  regularization  makes  the  method  significantly  more 
robust  to  signal  decorrelation.  However,  it  is  still  prone  to  decor¬ 
relation  within  each  window  especially  for  large  strain  rates.  In 
a  recent  work  [39],  a  displacement  field  is  first  calculated  by 
minimizing  phase  differences  in  correlation  windows  [21].  The 
strain  image  is  then  estimated  from  the  displacement  field  by  op¬ 
timizing  a  regularized  cost  function.  The  regularization  assures 
smooth  strain  image  calculation  from  the  noisy  displacement 
estimates. 

Dynamic  programming  (DP)  can  be  used  to  speed  the  op¬ 
timization  procedure  [40],  but  it  only  gives  integer  displace¬ 
ments.  Subsample  displacement  estimation  is  possible  [40],  but 
it  is  computationally  expensive,  particularly  if  subsample  accu¬ 
racy  is  needed  in  both  axial  and  lateral  directions.  Therefore, 
only  axial  subsample  displacement  is  calculated  [40].  In  ad¬ 
dition,  a  fixed  regularization  weight  is  applied  throughout  the 
image.  To  prevent  regions  with  high  local  decorrelation  from 
introducing  errors  in  the  displacement  estimation  one  should 
use  large  weights  for  the  regularization  term,  resulting  in  over¬ 
smoothing. 

In  this  paper,  we  present  two  novel  real-time  elastography 
methods  based  on  analytic  minimization  (AM)  of  cost  functions 
that  incorporate  similarity  of  echo  amplitudes  and  displacement 
continuity.  Similar  to  DP,  the  first  method  gives  subsample 
axial  and  integer  lateral  displacements.  The  second  method 
gives  subsample  2D  displacement  fields  and  2D  strain  fields. 
The  size  of  both  displacement  and  strain  fields  is  the  same  size 
as  the  RF-data  (i.e.,  the  methods  are  not  window  based  and  the 
displacement  and  strain  fields  are  calculated  for  all  individual 
samples  of  RF-data).  We  introduce  a  novel  regularization  term 
and  demonstrate  that  it  minimizes  displacement  underestima¬ 
tion  caused  by  smoothness  constraint.  We  also  introduce  the 
use  of  robust  statistics  implemented  via  iterated  reweighted 
least  squares  (IRLS)  to  treat  uncorrelated  ultrasound  data  as 
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outliers.  Finally,  for  the  first  time  to  the  best  of  our  knowledge 
we  introduce  the  use  of  Kalman  filtering  [41]  for  calculating 
strain  image  from  the  displacement  field.  Simulation  and  ex¬ 
perimental  results  are  provided  for  quantitative  validation.  The 
paper  concludes  with  a  clinical  pilot  study  utilizing  this  system 
for  monitoring  thermal  ablation  in  patients  with  liver  tumors. 

II.  Methods 

Assume  that  the  tissue  undergoes  a  deformation  and  let  I\ 
and  I2  be  two  images  acquired  from  the  tissue  before  and  after 
the  deformation.  Letting  Ii  and  I2  be  of  size  m  x  n  (Fig.  1), 
our  goal  is  to  find  two  matrices  A  and  L  where  the  (i,j) th  com¬ 
ponent  of  A(a-hj)  and  L(kj)  are  the  axial  and  lateral  motion 
of  the  pixel  (i,j)  of  Ii  (we  are  not  calculating  the  out-of-plane 
motion).  The  axial  and  lateral  strains  are  easily  calculated  by 
spatially  differentiating  A  in  the  axial  direction  (resulting  in  Aa) 
and  L  in  the  lateral  direction  (resulting  in  L\).  The  shear  strains 
(not  calculated  in  this  work)  can  also  be  easily  calculated  by 
spatially  differentiating  A  in  the  lateral  direction  (resulting  in 
Ai)  or  Lf  in  the  axial  direction  (resulting  in  La). 

In  this  section,  we  first  give  a  brief  overview  of  a  previous 
work  (DP)  that  calculates  integer  values  for  A  and  L.  We  then 
propose  ID  analytic  minimization  (AM)  as  a  method  that  takes 
the  integer  displacement  field  from  DP  and  refines  the  axial 
displacement  component.  We  then  introduce  2D  analytic  min¬ 
imization  (AM)  that  takes  the  integer  displacement  of  a  single 
RF-line  from  DP  and  gives  the  subsample  axial  and  lateral  dis¬ 
placement  fields  for  the  entire  image.  We  conclude  this  section 
by  presenting  a  technique  for  calculating  smooth  strain  field 
from  the  displacement  field  using  Kalman  filtering. 

A.  Dynamic  Programming  (DP) 

In  order  to  present  the  general  DP  formulation,  we  consider  a 
single  column  j  (an  RF-line)  in  Ii  (the  image  before  deforma¬ 
tion)  in  Fig.  1.  Let  m  and  n  be  the  length  of  the  RF-lines  and 
the  number  of  RF-lines  in  the  images  (Fig.  1).  Let  and  lt  de¬ 
note  the  axial  and  lateral  displacements  of  the  ith  sample  of  the 
RF-line  in  column  j.  In  DP  elastography  [40],  a  regularized  cost 
function  is  generated  by  adding  the  prior  of  displacement  conti¬ 
nuity  (the  regularization  term)  to  an  amplitude  similarity  term. 
The  displacement  continuity  term  for  column  j  is 

Rj  (  Q'i  1  Q'i—l  1  li—  l)  —  &i—  l)  A  li—  l)  (1) 

which  forces  the  displacements  of  the  sample  i  (i.e.,  Ui  and  k)  be 
similar  to  the  displacements  of  the  previous  sample  i  —  1  (i.e., 
at-i  and  k-i).  aa  and  ai  are  axial  and  lateral  regularization 
weights  respectively.  We  write  Rj(a^  1 1,  ai-i,  h-i)  to  indicate 
the  dependency  of  aL  and  ll  on  j.  The  regularized  cost  function 
for  column  j  is  then  generated  as  following: 

Cj  (O'li  ih  i) 

=  [h(ij)  ~h +  Ah )? 

.  f  Cj ( da ,  di^i  l)  -\-  Cj—i(da ,  di ,  i) 

+  nun  <  — - - - 

d^di  1  2 


where  da  and  di  are  temporary  displacements  in  the  axial  and 
lateral  directions  that  are  varied  to  minimize  the  term  in  the 
bracket.  After  calculating  Cj  for  i  =  2  . . .  m,  Cj  is  minimized 
at  i  =  rn  giving  am  and  lrn .  The  at  and  h  values  that  have 
minimized  the  cost  function  at  i  —  rn  are  then  traced  back  to 
i  —  1,  giving  integer  at  and  lt  for  all  samples  of  j  th  line.  The 
process  is  performed  for  the  next  line  j  +  1  until  the  displace¬ 
ment  of  the  whole  image  is  calculated.  The  2D  DP  method  gives 
integer  axial  and  lateral  displacement  maps.  In  [40],  we  per¬ 
formed  hierarchical  search  to  obtain  subsample  axial  displace¬ 
ment  (the  lateral  displacement  was  not  refined  to  subsample). 
DP  is  an  efficient  method  for  global  optimization  and  has  been 
used  extensively  in  many  applications  in  computer  vision  in¬ 
cluding  solving  for  optimal  deformable  models  [42].  In  the  next 
section,  we  propose  an  alternative  method  for  calculating  sub¬ 
sample  axial  displacement  which  is  both  faster  and  more  robust 
than  hierarchical  DP. 

B.  ID  Analytic  Minimization  (AM) 

Tissue  deformations  in  ultrasound  elastography  are  usually 
very  small  and  therefore  a  subsample  displacement  estimation  is 
required.  We  now  develop  a  method  that  analytically  minimizes 
a  regularized  cost  function  and  gives  the  refined  displacement 
field  following  the  work  presented  in  [16].  We  first  consider  a 
specialization  of  (2)  in  which  we  only  consider  refining  axial 
displacements  to  subsample  level. 

Having  the  integer  displacements  ai  and  lL  from  DP,  it  is  de¬ 
sired  to  find  A  a  ;  values  such  that  at  +  A  a  ;  gives  the  value  of 
the  displacement  at  the  sample  i  for  i  =  1 ...  m  ft,  at  and  A  ai 
correspond  to  line  j.  Hereafter,  wherever  the  displacements  cor¬ 
respond  to  the  j  th  line,  j  is  omitted  to  prevent  notation  clutter). 
Such  A  ai  values  will  minimize  the  following  regularized  cost 
function: 

Cj(Aai, . . . ,  A(2m) 

rn 

—  (hj)  ~  ^2(2 + a't + a^,  j + ii)]2 

i= 1 

+  H-  A &i  O'l—l  A(2^_x) 

A  ^i(^i  A  A cii  &itj-i  A(iy'_i)  ]}  (3) 

where  aa  >  0  and  ai  >  0  are  tunable  axial  and  lateral  regular¬ 
ization  weights  and  subscript  j  —  1  refers  to  the  previous  RF-line 
(adjacent  RF-line  in  the  lateral  direction). 

Substituting  I2(i  A  dt  +  Adi)  with  its  first-order  Taylor  ex¬ 
pansion  approximation  around  di ,  we  have 

m 

Cj'(Aai,  •  •  • ,  A(2m)  = 

i= 1 

—  12(1  A  cbi,j  A  h)  —  A  aij  +  /i)Aai)]2 
+  A  A  (hi  Q"i—1  A(2i_i) 

+  <y.i(ai  +  Aai  —  cbij-i  —  Aaij-i)2]}  (4) 


A  Rj  (ct"i  1 1i  1  da ,  dt) 


(2) 


where  H2  is  the  derivative  of  the  in  the  axial  direction.  The 
optimal  Aai  values  occur  when  the  partial  derivative  of  Cj  with 
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respect  to  Aa-b  is  zero.  Setting  ( dCj)/(dAai )  =  0  for  i  = 
1 ...  m  we  have 


Aa^  —  ^2e 


+  c^i)1 

a; 

+  ai&j-i, 

(5) 

-1 

0 

0  ••• 

0" 

2 

-1 

0  ••• 

0 

-1 

2 

_1  ... 

0 

(6) 

0 

0  -1 

1_ 

where  V2  =  diag^l  +  dij  +  /;) . .  ./^(m  +  dm,  j  +  k)), 

Aay  —  [Aaij  . . .  A(2m?ji']  ,  e  —  [ei . . .  cm]  ,  &i  —  I\  (i,  j ) 

I2 (i  +  dt,j  +  /t),  =  [ai  j  . . .  am ;j]T,  I  is  the  identity  matrix 
and  ay_i  is  the  total  displacement  of  the  previous  line  (i.e.,  when 
the  displacement  of  the  j  —  1th  line  was  being  calculated,  ay_i 
was  updated  with  a^-i  +  Aa^-i).  T>,  D  and  I  are  matrices  of 
size  m  x  m  and  Aa,  e  and  a  are  vectors  of  size  m. 

Comparing  ID  AM  [as  formulated  in  (5)]  and  2D  DP,  they 
both  optimize  the  same  cost  function.  Therefore,  they  give  the 
same  displacement  fields  (up  to  the  refinement  level  of  the  DP). 
In  the  next  two  subsections,  we  will  further  improve  ID  AM. 

1)  Biasing  the  Regularization :  The  regularization  term 
aa(a,i  +  A Oi  —  at-i  —  Aa;_i)2  penalizes  the  difference 
between  at  +  Aat  and  at-i  +  Aat_i,  and  therefore  can 
result  in  underestimation  of  the  displacement  field.  Such 
underestimation  can  be  prevented  by  biasing  the  regulariza¬ 
tion  by  e  to  aa(a,i  +  A —  at-i  —  Aa.t_i  —  e)2,  where 
e  =  (am  —  ai)/(m  —  1)  is  the  average  displacement  difference 
(i.e.,  average  strain)  between  samples  i  and  i  —  1.  An  accurate 
enough  estimate  of  dm  —  d\  is  known  from  the  previous 
line.  With  the  bias  term,  the  right-hand  side  of  (5)  becomes 
T>e  —  (aaD  +  &iI)aLj  +  ai(aLj-i  +  Aa^-i)  +  b  where  the 
bias  term  is  b  =  aa[— t  0  ...  0  e]T  (only  the  first  and  the  last 
terms  are  nonzero)  and  all  other  terms  are  as  before.  In  the  other 
words,  except  for  the  first  and  the  last  equations  in  this  system, 
all  other  m  —  2  equations  are  same  as  (5). 

Equation  (5)  can  be  solved  for  Aa^  in  4m  operations  since 
the  coefficient  matrix  I2  +  aaD  +  ail  is  tridiagonal.  Utilizing 
its  symmetry,  the  number  of  operations  can  be  reduced  to  2m. 
The  number  of  operations  required  for  solving  a  system  with 
a  full  coefficient  matrix  is  more  than  m3/ 3,  significantly  more 
than  2m. 

2)  Making  Elastography  Resistant  to  Outliers:  Even  with 
pure  axial  compression,  some  regions  of  the  image  may  move 
out  of  the  imaging  plane  and  decrease  the  decorrelation.  In  such 
parts  the  weight  of  the  data  term  in  the  cost  function  should  be 
reduced.  The  data  from  these  parts  can  be  regarded  as  outliers 
and  therefore  a  robust  estimation  technique  can  limit  their  effect. 
Before  deriving  a  robust  estimator  for  Ad,  we  rewrite  (4)  as 

m 

C(Ad)  =  ^(rO  +  ii(Ad)  (7) 

i= 1 

where  rt  —  Ii(i)  —  /2O'  +  di)  —  I2(i  +  di)Adi  is  the  residual, 
p(ri)  =  r\  and  R  is  the  regularization  term.  The  M-estimate  of 
Ad  is  Ad  =  argminAd{^Hi  p(n)  +  A!(Ad)}  where  p(n) 


is  a  robust  loss  function  [43].  The  minimization  is  solved  by 
setting  (OC)/(dAdi)  —  0 


An) 


dr 

dA  di 


+ 


dR(Ad) 

dAdi 


=  0. 


(8) 


A  common  next  step  [44]  is  to  introduce  a  weight  function  w, 
where  w{ri).ri  —  pf{ri).  This  leads  to  a  process  known  as  “iter¬ 
atively  reweighted  least  squares”  (IRLS)  [45],  which  alternates 
steps  of  calculating  weights  w(ri)  for  =  1 . . .  m  using  the 
current  estimate  of  Ad  and  solving  (8)  to  estimate  a  new  Ad 
with  the  weights  fixed.  Among  many  proposed  shapes  for  w(  • ), 
we  compared  the  performance  of  Huber  [44],  [43] 

,  ,  f  1  Ird  <T 

W{ri)  =  1  R  W  >  T  <9> 


and  Cauchy  [45] 

w(n)  =  1 + <kfr-  (10) 

functions  and  discovered  that  the  more  strict  Cauchy  function 
(which  decreases  with  inverse  of  the  square  of  the  residual)  is 
more  suitable  in  our  application.  To  better  discriminate  outliers, 
we  calculate  the  residuals  r  t  at  linear  interpolation  of  the  integer 
sample  displacements  provided  by  DP.  With  the  addition  of  the 
weight  function,  (8)  becomes 

(wT22  +  aaT>  +  ttJ!)  Aa j  =  wl'2e 

“(ftaD  +  +  ai&j-i  +  b  (11) 

where  w  =  diag(u?(ri) . . .  w(rm)).  This  equation  will  con¬ 
verge  to  a  unique  local  minimum  after  few  iterations  [45].  The 
convergence  speed  however  depends  on  the  choice  of  T,  which 
in  this  work  is  defined  manually.  Since  the  Taylor  approxima¬ 
tion  gives  a  local  quadratic  approximation  of  the  original  non¬ 
quadratic  cost  function,  the  effect  of  higher  orders  terms  in¬ 
crease  if  A dj  is  large.  Assuming  that  DP  gives  the  correct  dis¬ 
placements,  1 1  Acij  1 1  oc  <  e  where  ||  •  is  the  infinity  norm 
and  e  <  0.5.  In  practice,  however,  e  <  0.5  because  the  linear 
interpolation  of  the  DP  displacements  (which  is  very  close  to 
the  correct  displacement)  is  used  to  calculate  the  residuals  rt. 
Therefore,  a  small  value  can  be  assigned  to  T  in  ID  AM  pro¬ 
vided  that  DP  results  are  trusted. 

The  coefficient  matrix  Q  —  wl2  +  aa  D  +  ail  in  (11)  is  the 
Hessian  of  the  cost  function  C  whose  minimum  is  sought.  This 
matrix  is  strictly  diagonally  dominant  (i.e.,  \qvl\  >  l^-l 

for  all  i  where  qtj  is  the  i,j th  element  of  Q),  symmetric  and  all 
diagonal  entries  are  positive.  Therefore,  it  is  positive  definite, 
which  means  that  setting  the  gradient  of  C  to  zero  results  in  the 
global  minimum  of  C  (not  in  a  saddle  point,  a  local  maximum 
or  a  local  minimum).  All  of  the  ID  AM  results  presented  in  this 
work  are  obtained  with  one  iteration  of  the  above  equation. 

ID  AM  takes  the  integer  axial  and  lateral  displacement  fields 
from  DP  and  gives  refined  axial  displacement.  It  inherits  the 
robustness  of  DP  and  adds  more  robustness  when  calculating  the 
fine  axial  displacements  via  IRLS.  However,  there  are  redundant 
calculations  in  this  method  which  are  eliminated  in  2D  AM  as 
described  next. 
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C.  2D  Analytic  Minimization  (AM) 

In  2D  AM,  we  modify  (2)  to  calculate  subsample  axial  and 
lateral  displacement  fields  simultaneously.  The  outline  of  our 
proposed  algorithm  is  as  follows. 

1)  Calculate  the  integer  axial  and  lateral  displacements  of  one 
or  more  seed  RF-lines  (preferably  in  the  middle  of  the 
image)  using  DP  [(2)] .  Calculate  the  linear  interpolation  of 
the  integer  displacements  as  an  initial  subsample  estimate. 

2)  Calculate  subsample  axial  and  lateral  displacements  of  the 
seed  RF-line  using  2D  AM,  as  explained  below.  Add  the 
subsample  axial  and  lateral  displacements  to  the  initial  es¬ 
timate  to  get  the  displacement  of  the  seed  line. 

3)  Propagate  the  solution  to  the  right  and  left  of  the  seed 
RF-line  using  the  2D  AM  method,  taking  the  displacement 
of  the  previous  line  as  the  initial  displacement  estimate. 

Benefits  of  2D  AM  are  two-fold.  First  it  computes  subsample 
displacements  in  both  axial  and  lateral  directions.  Lateral  strain 
contains  important  information  from  tissue  structure  that  is  not 
available  from  axial  strain  [31],  [46],  [47].  Second,  it  is  only 
required  to  calculate  the  displacement  of  a  single  line  using  DP 
(the  seed),  eliminating  the  need  to  have  the  integer  displacement 
map  for  the  entire  image.  This  is  significant  as  in  the  ID  AM 
method,  the  initial  step  to  calculate  the  2D  integer  displacements 
using  DP  takes  about  10  times  more  than  the  ID  AM. 

Assume  that  initial  displacement  estimates  in  the  axial  di¬ 
rection,  ai,  and  in  the  lateral  direction,  lt,  are  known  for  all 
i  —  1 . . .  m  samples  of  an  RF-line.  Note  that  at  and  lt  are  not 
integer;  for  the  seed  line  they  are  the  linear  interpolation  of  the 
integer  DP  displacements  and  for  the  rest  of  the  lines  are  the 
displacement  of  the  previous  line.  It  is  desired  to  find  A  a;  and 
Ali  values  such  that  the  duple  (at  +  A ll  +  A Ip  gives  the 
axial  and  lateral  displacements  at  the  sample  i.  Such  (A  di,  A  ap 
values  will  minimize  the  following  regularized  cost  function: 

Cy(Aai, . . . ,  Aam,  AZi, . . . ,  A Zm) 

m 

=  _  ^2(2 + ai + j + ii + Aii)]2 

i= 1 

+  a(ai  +  Aai  -  ai- 1  -  Aai-i )2 

+  Pa(li  +  AZ i  ~  h- 1  —  Ali-i)2 

+  Pi(k  +  Alt  —  h,j-i)2}  (12) 

where  is  the  ith  sample  on  the  j th  RF-line.  Since  we 

perform  the  calculations  for  one  RF-line  at  a  time,  we  dropped 
the  index  j  to  simplify  the  notations:  at,  /t,  Aat,  and  Alt  are 
ahj ,  lhj ,  Aahj ,  and  A lhj.  lij-i  is  the  lateral  displacement  of 
the  previous  RF-line  (note  that  hj-i  is  the  total  lateral  dis¬ 
placement  of  the  previous  line,  i.e.,  when  the  displacement  of 
the  j  —  1th  line  was  being  calculated,  kj-i  was  updated  with 
hj-i  +  AZij^i).  Since  in  the  first  iteration  at  and  k  (the  ini¬ 
tial  displacement  estimates)  are  in  fact  the  displacements  of  the 
previous  RF-line,  for  the  first  iteration  we  have  hj-i  —  h.  This 
simplifies  the  last  term  in  the  right-hand  side  to  (3[A l2.  The  reg¬ 
ularization  terms  are  a,  j3a  and  (3[\  a  determines  how  close  the 
axial  displacement  of  each  sample  should  be  to  its  neighbor  on 
the  top  and  [3a  and  [3[  determine  how  close  lateral  displacement 
of  each  sample  should  be  to  its  neighbors  on  the  top  and  left 


(or  right  if  propagating  to  the  left).  If  the  displacement  of  the 
previous  line  is  not  accurate,  it  will  affect  the  displacement  of 
the  next  line  through  the  last  term  in  the  right-hand  side  of  (12). 
Although  its  effect  will  decrease  exponentially  with  j,  it  will 
propagate  for  few  RF-lines.  Therefore,  we  set 


tf  = 


Pi 

1  + 


(13) 


to  prevent  such  propagation  where  rij-i  is  the  residual  asso¬ 
ciated  with  the  displacement  of  the  Ah  sample  of  the  previous 
line.  A  large  residual  indicates  that  the  displacement  is  not  accu¬ 
rate  and  therefore  its  influence  on  the  next  line  should  be  small, 
which  is  realized  via  the  small  weight  (3[.  This  is,  in  principle, 
similar  to  guiding  the  displacement  estimation  based  on  a  data 
quality  indicator  [48].  The  effect  of  the  tunable  parameters  a, 
[3 a  and  (3i  is  studied  in  Section  III.  Writing  the  2D  Taylor  ex¬ 
pansion  of  the  data  term  in  (12)  around  (i  +  at,j  +  Ip 


I2P  +  &i  +  Aa^j  +  li  +  AZi) 

~  I2P  +  ai,j  +  li)  +  Aail^  a  +  Al-J^i  (14) 


where  V2  a  and  H2  t  are  the  derivatives  of  the  I2  at  point  (i  + 
ai,j  +  li)  in  the  axial  and  lateral  directions  respectively.  Note 
that  since  the  point  (i  +  at ,  j  +  li)  is  not  on  the  grid  (at  and 
li  are  not  integer),  interpolation  is  required  to  calculate  If2  a 
and  I f2  v  We  propose  a  method  in  Section  II-C1  that  eliminates 
the  need  for  interpolation.  The  optimal  (A Alp  values  occur 
when  the  partial  derivatives  of  Cj  with  respect  to  both  Aat  and 
All  are  zero.  Setting  (dCj)/(dAap  =  0  and  (dCj)/(dAlp  = 
0  for  i  —  1 . . .  m  and  stacking  the  2 m  unknowns  in  Ad  = 
[Aai  Zi  Aa2  A/2  . . .  Aam  A lm]T  and  the  2m  initial  estimates 


L  _L  _L  ^  ^ 

in  d  =  [ai  Zi  a2  I2  •  • 
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J 

T  we  have 
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where  V2  —  diag(0,  (3[,  0,  (3[, . . . ,  0,  (3[)  is  a  diagonal  matrix  of 
size  2m  x  2m,  Xf2  =  diag(3/2(l) . . .  3/2(m))  is  a  symmetric 
tridiagonal  matrix  of  size  2m  x  2m  with 


= 


v  1 

12,a 


Tf  Tf 

12,a12,l 

Tf  Jf  Tf  2 

\-I2,aI2,l  L2 ,1 


(16) 


blocks  on  its  diagonal  entries  where  I2  a  and  I2  t  are  the  deriva¬ 
tives  of  the  I 2  at  point  (i  +  at,j  +  Ip  in  the  axial  and  lateral 
directions 


1'2  =  diag(/'ia(l),  7^(1),  1^(2),  I'n{2) . . . 

(17) 
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where  I2  a(i)  and  are  calculated  at  point  (i  +  at,j  +  li), 

and  e  =  [ei  ei  e2  e2  . . .  em]T,  ei  =  h(i,j)  - h(i  +  ai,j  +  k). 

We  make  four  modifications  to  (15).  First,  we  take  into  ac¬ 
count  the  attenuation  of  the  ultrasound  signal  with  depth.  As  the 
signal  gets  weaker  with  depth,  the  first  term  in  the  right-hand 
side  of  (15)  (I2 ft)  gets  smaller.  This  results  in  increasing  the 
share  of  the  regularization  term  in  the  cost  Cj  and  therefore 
over-smoothing  the  bottom  of  the  image.  The  attenuation  of 
the  ultrasound  signal  [49]  reflected  from  the  depth  d  is  ((d)  = 
e-2iog(io)at/0d/20  where  af  [ s  the  frequency  dependent  attenu¬ 
ation  coefficient  of  tissue  and  is  equal  to  0.63  dB/cm/MHz  for 
fat  [49],  /o  is  the  center  frequency  of  the  wave  (in  MHz)  and  d 
is  in  cm.  Having  the  exponential  attenuation  equation,  the  atten¬ 
uation  level  at  sample  i  will  be 

154Q  X  lQ2at/Q  log(10) 

(i  =  X~\  X  —  e  20/s:Sc;106  5  i  —  1  .  .  .  m  (18) 

where  1540  x  102  is  the  speed  of  sound  in  tissue  (in  cm/sec) 
and  fs  is  the  sampling  rate  of  the  ultrasound  system  (in  MHz). 
This  is  assuming  that  the  TGC  (time  gain  control)  is  turned 
off.  Otherwise,  the  TGC  values  should  be  taken  into  account 
in  this  equation.  Let  the  2m  x  2m  diagonal  matrix  Z  be  Z  = 
diag(£i,  Ci,  C2,  (2  •  •  •  Cm,  Cm)-  To  compensate  for  the  attenua¬ 
tion,  we  multiply  the  T>\  and  V2  matrices  in  (15)  by  Z,  and 
therefore  reduce  the  regularization  weight  with  depth.  As  we 
will  show  in  Sections  III  and  IV,  the  regularization  weight  can 
vary  substantially  with  no  performance  degradation.  Therefore 
approximate  values  of  the  speed  of  sound  and  attenuation  co¬ 
efficient  will  suffice.  Second,  we  add  a  bias  term  in  the  regu¬ 
larization  similar  to  the  ID  case.  Here  we  only  bias  the  axial 
displacement  since  the  difference  between  the  lateral  displace¬ 
ments  of  the  points  on  a  RF-line  is  very  small,  usually  less  than 
4  RF-lines.  Third,  we  exploit  the  fact  that,  because  the  tissue 
is  in  contact  with  the  ultrasound  probe,  the  axial  displacement 
of  the  top  of  the  image  is  zero  relative  to  the  probe  (the  lateral 
displacement  of  the  top  of  the  image  is  not  zero  as  tissue  might 
slip  under  the  probe).  Therefore,  we  enforce  the  axial  displace¬ 
ment  of  the  first  sample  to  be  zero  by  changing  the  first  row  of 
V 1,  Z2,  and  I2  .  Fourth,  we  make  the  displacement  estimation 
robust  via  IRLS  using  the  Cauchy  function  (10).  Similar  to  ID 
AM,  T  is  selected  manually.  For  the  first  (seed)  RF-line,  a  small 
value  can  be  selected  for  T  if  DP  results  are  trusted.  For  the  next 
lines,  the  value  of  Ad  determines  the  accuracy  of  the  Taylor  ex¬ 
pansion  14:  for  a  small  Ad,  the  residuals  of  the  inliers  are  small 
and  therefore  a  small  T  can  be  chosen,  while  for  a  large  Ad  the 
inliers  might  give  large  residuals  and  therefore  a  large  value  for 
T  is  required.  Since  the  tissue  motion  is  mostly  continuous,  Ad 
mostly  depends  on  the  lateral  sampling  of  the  image  (i.e.,  the 
number  of  A-line  per  cm).  Therefore  if  many  A-lines  are  given 
per  cm  of  the  image  width,  a  small  value  of  T  will  give  the  op¬ 
timum  results.  Since  the  amplitude  of  signal  is  decreasing  due 
to  attenuation,  we  decrease  the  IRLS  parameter  T  with  depth 
by  multiplying  it  with  Q  at  each  sample  i.  With  these  modifica¬ 
tions,  (15)  becomes 

(' Wlr2 2  +  ZV1  +  ZX>2)Ad  =  Wl/e  -  ZDlCl  +  s  (19) 

where  W  =  diag(0,u>(ri),u>(?-2),i(;(7-2) . .  -w(rm),w(rm)) 
(i.e.,  W2;,2i  =  W2t— i;2i— i  =  w(n)  for  i  =  1 . . .  m  except  for 
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Wi,i  =  0  which  guarantees  the  displacement  of  the  first  sample 
to  be  zero)  is  the  weight  function  determined  by  the  residuals 
n  =  h(i,j)  -  [h(i  +  di,j  +  Oi)  +  Al dil'2  z  +  A w  is 
as  before  (10),  the  bias  term  s  is  a  vector  of  length  2m  whose 
all  elements  are  zero  except  the  2m  —  1th  element: S2m-i  —  ae, 
and  e  =  (drn  —  di)/(m  —  1)  is  as  before.  Similar  to  (1 1),  the  co¬ 
efficient  matrix  Q  —  WZ2  +  ZX>i  +  ZX>2  is  strictly  diagonally 
dominant,  symmetric  and  all  the  diagonal  entries  are  positive. 
Therefore  Q  is  positive  definite  which  means  that  solving  (19) 
results  in  the  global  minimum  of  the  cost  function  C.  The  up¬ 
dated  displacement  field  (axial  and  lateral)  will  be  d  +  Ad. 

Equation  (19)  can  be  solved  for  Ad  in  9  m  operations  since 
the  coefficient  matrix  WZ2  +  ZPi  +  ZZ>2  is  pentadiagonal 
and  symmetric.  This  number  is  again  significantly  less  than 
((2m)3)/(3),  the  number  of  operations  required  to  solve  a  full 
system. 

1 )  Inverse  Gradient  Estimation:  With  the  subsample  initial 
displacement  field,  the  Taylor  expansion  should  be  written 
around  off-grid  points,  which  requires  calculation  of  image 
gradient  at  these  points  [matrices  Z2  and  Z2  in  (19)].  In 
Fig.  2(a),  this  is  equivalent  to  calculating  gradient  of  I 2  on  the 
off-grid  marks.  There  are  two  disadvantages  associated  with 
this:  1)  it  requires  interpolation  of  the  gradients,  and  2)  the 
image  gradient  should  be  recalculated  after  each  iteration.  As 
proposed  by  [44],  [50],  image  gradient  can  be  instead  calculated 
at  on-grid  locations  on  image  1  in  the  following  way. 

Consider  two  problems:  1)  to  find  the  matches  for  grid  points 
on  Ji  having  the  initial  off-grid  estimates  on  /2,  and  2)  to  find 
the  matches  for  the  off-grid  points  on  /2  having  the  initial  grid 
estimates  on  /1 .  For  both  problems,  I 2  values  must  be  interpo¬ 
lated  on  the  off-grid  values.  However,  the  second  problem  does 
not  require  interpolation  of  the  image  gradient  since  the  Taylor 
expansion  is  written  around  grid  points  of  Ii  [Fig.  2(b)].  It  is 
shown  in  [51]  that  the  two  techniques  converge  to  the  same  re¬ 
sults.  Therefore,  on  one  hand  inverse  gradient  calculation  is  both 
faster  and  easier  to  implement,  and  on  the  other  hand  it  causes 
no  performance  degradation.  Exploiting  this,  (19)  becomes 

(WZ[2  +  ZDi  +  ZP2)  Ad  =  Wl/e  -  ZPxd  +  s  (20) 

where  T2  and  Z(  are  now  calculated  on  the  grid  points  of  image 
1. 

All  the  2D  AM  results  presented  in  this  work  are  obtained 
using  (20).  For  the  seed  line  where  the  initial  estimate  might 
be  inaccurate,  this  equation  is  iterated  multiple  times  (about  10 
times).  For  all  other  lines,  this  equation  is  iterated  only  once. 

D.  Strain  Estimation  Using  Kalman  Filter 

Strain  estimation  requires  spatial  derivation  of  the  displace¬ 
ment  field.  Since  differentiation  amplifies  the  signal  noise,  least 
squares  regression  techniques  are  commonly  used  to  obtain  the 
strain  field.  Adjacent  RF-lines  are  usually  processed  indepen¬ 
dently  in  strain  calculation.  However,  the  strain  value  of  each 
pixel  is  not  independent  from  the  strain  value  of  its  neighboring 
pixels.  The  only  exception  is  the  boundary  of  two  tissue  types 
with  different  mechanical  properties  where  the  strain  field  is  dis¬ 
continuous.  We  use  the  prior  of  piecewise  strain  continuity  via 
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(a)  (W 

Fig.  2.  (a)  In  I2  the  initial  estimates  (in  black)  are  updated  by  the  arrows  (three  components  of  Ad)  to  new  estimates  (in  red)  after  an  iteration  of  2D  AM.  To 
find  Ad  using  (19),  it  is  required  to  calculate  image  gradient  at  the  off-grid  initial  estimate  locations  (in  black)  on  I2.  (b)  Schematic  plot  of  two  RF-data  A  and 
I2 ,  each  sampled  at  three  locations  (black  dots).  The  black  dashed-dotted  arrow  shows  A  a  of  the  sample  on  I±  (ignoring  the  regularization  term)  which  requires 
calculating  the  gradient  on  I2  at  an  off-grid  location.  The  blue  dashed  arrow  shows  A  a  of  an  off-grid  sample  on  I2  (ignoring  the  regularization  term)  which  requires 
calculating  the  gradient  on  A  at  an  on-grid  location.  Ignoring  second-order  derivatives,  the  length  of  the  two  arrows  is  equal,  (a)  Three  samples  on  A  (left)  and 
corresponding  matches  on  I2  (right),  (b)  Inverse  gradient  estimation. 


a  Kalman  filter  to  improve  the  quality  of  strain  estimation.  Al¬ 
though  locations  with  strain  discontinuity  are  limited,  we  will 
develop  a  technique  to  take  such  discontinuities  into  account. 

We  first  calculate  the  strain  using  least  squares  regression. 
Each  RF-line  is  first  differentiated  independently:  for  each 
sample  i ,  a  line  is  fitted  to  the  displacement  estimates  in  a 
window  of  length  2k  +  1  around  i ,  i.e.,  to  the  samples  i  —  k 
to  i  +  k.  The  slope  of  the  line,  Zij,  is  calculated  as  the  strain 
measurement  at  i.  The  center  of  the  window  is  then  moved  to 
i  +  1  and  the  strain  value  a+i  j  is  calculated.  We  reuse  over¬ 
lapping  terms  in  calculation  of  Zij  and  Zi+ ij,  and  therefore 
the  running  time  is  independent  of  the  window  length  2k  +  1. 
Having  z%j  for  i  —  1 . . .  m  and  j  —  1 ...  n,  we  propose  the 
following  algorithm  based  on  Kalman  filter  to  take  into  account 
the  prior  of  strain  continuity. 

Zij  are  the  noisy  measurements  of  the  underlying  strain  field 
tij.  Since  the  ztj  values  are  calculated  using  axial  windows, 
we  apply  the  Kalman  filter  in  the  lateral  direction.  Let  ry  be  the 
Gaussian  process  noise  and  Sj  be  the  Gaussian  measurement 
noise  to  be  removed.  We  have  [52],  [41] 

tij  —  A  j-i  +  nj  (21) 

Z-ij  —  tij  +  S'ij-  (22) 

Let  e~j  (note  the  super  minus)  be  our  a  priori  strain  estimate 
from  the  process  prior  to  step  j  [i.e.,  from  the  (21)]  and  tij  be 
our  a  posteriori  strain  estimate  at  step  j  given  measurement  Zj . 
Let  also  the  variances  of  t~.-  and  be  respectively  p~  and  p. 
The  time  update  (i.e.,  prior  estimation)  equations  will  be  [41] 


kj  =  ki- 1  (23) 

Pi,j  =  Pi, 3-1  +  ar  (24) 


where  a 2  is  the  variance  of  the  process  noise  r.  Vi,j-i  is  initial- 
ized  to  zero  for  the  first  sample  j  —  1 .  The  measurement  update 
equations  will  be  [41] 


id  -  hj  +  -  ’  2  ~ 

Pij  + 


Pi, 3  =  1 


Pi, 3 


Pi, 3  +  a 


Pi 


2  /  ^ 


(25) 

(26) 


where  is  the  variance  of  the  measurement  noise  5.  Note  that 
since  both  the  state  tij  and  measurement  zij  are  scalars,  all  the 
update  equations  only  require  scalar  operations.  We  estimate 
and  a'l  as  following.  Let  the  mean  (calculated  using  a  Gaussian 
kernel  of  standard  deviation  of  cfg  —  0.6  sample)  of  the  strain 
values  in  3  x  3  blocks  around  samples  (i,j  —  1)  and  (i,j)  be 
fjj-i  and  pj,  respectively.  Then  is  [52] 


j  =  (fij- 1  -  Hj)2.  (27) 


This  is  a  reasonable  estimate  of  a %  as  it  tries  to  capture  the  dif¬ 
ference  between  pixel  values  at  adjacent  RF-lines.  If  the  dif¬ 
ference  between  the  mean  strain  values  is  high,  less  weight  is 
given  to  the  a  priori  estimate.  This  space-variant  estimation  of 
the  model  noise  provides  a  better  match  to  local  variations  in  the 
underlying  tissue  leading  to  a  greater  noise  reduction,  a f  is  the 
variance  of  zL_j  measurements  in  the  entire  image  and  is  con¬ 
stant  throughout  the  image. 

The  strain  estimation  algorithm  can  be  summarized  as 
following. 

1)  Perform  least  squares  regression  in  the  axial  direction  for 
each  RF-line.  Generate  a  (noisy)  strain  image  Z  whose 
pixel  i,j  is  Zij.  This  step  ensures  continuity  in  the  axial 
direction. 

2)  Apply  the  Kalman  filter  to  the  noisy  strain  image  Z  in  the 
lateral  direction.  Generate  a  (denoised)  strain  image  whose 


RIVAZ  et  al. :  REAL-TIME  REGULARIZED  ULTRASOUND  ELASTOGRAPHY 


935 


pixel  i,  jf  is  This  step  ensures  continuity  in  the  lateral 
direction. 

Both  steps  are  applied  once  and  are  not  iterated.  We  show  in  the 
experimental  results  how  the  Kalman  filter  removes  the  noise 
from  the  strain  image  with  minimal  blurring,  owing  to  the  model 
noise  update  (27). 


III.  Simulation  Results 

Field  II  [53]  and  ABAQUS  (Providence,  RI)  software  are 
used  for  ultrasound  simulation  and  for  finite  element  simula¬ 
tion.  Many  scatterers  are  distributed  in  a  volume  and  an  ultra¬ 
sound  image  is  created  by  convolving  all  scatterers  with  the 
point  spread  function  of  the  ultrasound  and  adding  the  results 
using  superposition.  The  phantom  is  then  meshed  and  com¬ 
pressed  using  finite  element  simulation,  giving  the  3D  displace¬ 
ment  of  each  node  of  the  mesh.  The  displacement  of  each  scat- 
terer  is  then  calculated  by  interpolating  the  displacement  of  its 
neighboring  nodes.  Scatterers  are  then  moved  accordingly  and 
the  second  ultrasound  image  is  generated.  The  displacement  and 
strain  fields  are  then  calculated  using  the  AM  methods  and  are 
compared  with  the  ground  truth.  The  unitless  metric  signal-to- 
noise  ratio  (SNR)  and  contrast  to  noise  ratio  (CNR)  are  also  cal¬ 
culated  to  assess  the  performance  of  the  AM  method  according 
to 


SNR  =  - 
a 


(28) 


where  st  and  si>  are  the  spatial  strain  average  of  the  target  and 
background,  af  and  of  are  the  spatial  strain  variance  of  the 
target  and  background,  and  s  and  a  are  the  spatial  average  and 
variance  of  a  window  in  the  strain  image,  respectively. 

The  parameters  of  the  ultrasound  probe  are  set  to  mimic  com¬ 
mercial  probes.  The  probe  frequency  is  7.27  MHz,  the  sampling 
rate  is  40  MHz  and  the  fractional  bandwidth  is  60%.  A  Hanning 
window  is  used  for  apodization,  the  single  transmit  focus  is  at 
22.5  mm,  equi-distance  receive  foci  are  from  5  to  45  mm  at  each 
5  mm,  the  transmit  is  sequential,  and  the  number  of  active  ele¬ 
ments  is  64. 

Two  simulated  phantoms  are  generated.  The  first  phantom 
is  50  x  10  x  55  mm  and  the  second  one  is  36  x  10  x  25  mm. 
Respectively  5  x  10°  and  1.4  x  10°  scatterers  with  Gaussian 
scattering  strengths  [54]  are  uniformly  distributed  in  the  first  and 
second  phantom,  ensuring  more  than  10  scatterers  [55]  exist  in 
a  resolution  cell. 

The  mechanical  properties  of  both  phantoms,  required  for  fi¬ 
nite  element  simulation,  is  assumed  to  be  isotropic  and  homoge¬ 
neous.  The  first  phantom  is  uniform  while  the  second  phantom 
contains  a  circular  hole  filled  with  blood  that  can  move  out-of- 
plane,  simulating  a  blood  vessel  in  tissue  [Fig.  7(a)].  The  scat¬ 
terers  are  distributed  in  the  vessel,  also  with  the  same  intensity 
and  distribution  as  the  surrounding  material.  A  uniform  com¬ 
pression  in  the  x  direction  is  applied  and  the  3D  displacement 
field  of  phantoms  is  calculated  using  ABAQUS.  The  Poisson’s 
ratio  is  set  to  v  —  0.49  in  both  phantoms  to  mimic  real  tissue 
[56],  [57]  which  causes  the  phantoms  to  deform  in  x  and  y  di¬ 
rections  as  a  result  of  the  compression  in  the  x  direction. 


The  first  phantom  undergoes  compressions  in  the  x  direction 
to  achieve  strain  levels  of  1%-10%.  Fig.  3  shows  the  SNR  of 
the  axial  strain  of  the  ID  AM  and  2D  AM  methods  [the  window 
for  SNR  calculation  covers  the  entire  strain  image  in  (a)  and  (f)]. 
The  sharp  drop  of  the  SNR  with  strain  in  graph  (a)  is  mainly  due 
to  the  strain  underestimation  in  the  bottom  part  of  the  image.  It 
can  be  explained  as  following.  The  unbiased  regularization  term 
tries  to  force  constant  displacement  [dashed-dotted  red  line  in 
(b)].  Assuming  an  ideal  noiseless  case  where  the  data  term  gives 
a  smooth  ramp  displacement  [dashed  black  line  in  (b)],  mini¬ 
mizing  the  cost  function  (which  is  the  summation  of  the  data  and 
the  regularization  terms)  will  underestimate  the  displacement  at 
the  two  ends  [solid  blue  line  in  (b)] .  This  underestimation  decays 
exponentially  moving  towards  the  center  of  the  image.  This  ar¬ 
tifact  is  shown  in  the  simulation  experiment  at  2%  and  6%  strain 
levels  in  (c)  and  (d).  Since  we  exploit  the  fact  that  the  axial  dis¬ 
placement  of  the  first  sample  is  zero  (Section  II-C),  the  under¬ 
estimation  does  not  happen  in  the  top  of  the  image.  Biasing  the 
regularization  prevents  this  artifact,  as  is  shown  in  (c)  and  (d). 
The  AM  method  with  or  without  the  bias  term  gives  the  same 
result  away  from  the  bottom  of  the  image:  part  (e)  shows  that  if 
we  ignore  300  (5.8  mm)  samples  at  the  bottom  of  the  image,  the 
SNR  will  not  drop  sharply  unlike  in  part  (a).  Part  (f)  shows  the 
SNR  of  the  AM  methods  with  biased  regularization  calculated 
in  the  entire  image.  The  SNR  at  1%  strain  in  parts  (e)  and  (f) 
is  the  same.  At  higher  strain  levels,  the  strain  underestimation 
propagates  more  into  the  middle  of  the  image,  and  therefore  the 
SNR  decreases  at  higher  strain  levels  in  graph  (e).  Part  (e)  shows 
2D  AM  gives  slightly  better  axial  strain  compared  to  ID  AM. 
IRLS  slightly  increases  the  SNR.  However,  we  will  see  in  the 
simulation  results  of  the  second  phantom  that  in  the  presence  of 
outliers  significant  improvement  in  SNR  and  CNR  is  achieved 
using  IRLS. 

The  SNR  of  the  lateral  strain  field  is  much  lower  than  that  of 
the  axial  strain  field  (Fig.  4).  Unbiased  regularization  gives  the 
lowest  SNR,  mainly  due  to  artifacts  in  the  bottom  of  the  image. 
Similar  to  the  axial  strain,  the  SNR  improves  as  300  samples 
from  the  bottom  of  image  are  omitted  from  the  SNR  calculation 
(results  not  shown). 

The  effect  of  the  regularization  weights  on  bias  and  variance 
of  the  axial  strain  image  at  2%  ground  truth  axial  strain  is  shown 
in  Fig.  5.  The  blue  curves  show  the  bias  and  variance  of  the  en¬ 
tire  strain  image  obtained  with  unbiased  regularization.  It  shows 
the  tradeoff  between  the  bias  and  variance:  increasing  the  regu¬ 
larization  weight  increases  the  bias  and  decreases  the  variance. 
The  variance  starts  to  increase  at  a  ^  12  which  is  caused  by  the 
underestimation  of  the  strain  at  the  bottom  of  the  image  [the  ar¬ 
tifact  in  Fig.  3(c)].  If  we  exclude  the  bottom  300  samples  of  the 
strain  image  from  the  bias  and  variance  calculation  (the  black 
dashed  curve),  we  see  a  consistent  drop  of  variance  as  a  is  in¬ 
creased.  The  black  curves  show  the  bias  and  variance  of  the  en¬ 
tire  strain  image  obtained  with  biased  regularization.  Biasing  the 
regularization  causes  the  bias  to  decrease  as  the  regularization 
weight  a  is  increased  which  is  a  nonstandard  behavior.  It  can  be 
explained  by  the  simple  ground  truth  strain  field  which  is  uni¬ 
form,  exactly  what  the  regularization  term  is  trying  to  achieve. 
Even  in  the  unbiased  case,  only  the  bias  of  the  bottom  part  of  the 
strain  field  increases  as  a  is  increased  (i.e.,  in  the  bias  plot,  the 


936 


IEEE  TRANSACTIONS  ON  MEDICAL  IMAGING,  VOL.  30,  NO.  4,  APRIL  2011 


(d) 


<e> 


(0 


Fig.  3.  Axial  strain  estimation  in  the  first  simulated  phantom,  (a)  The  SNR  values  corresponding  to  the  unbiased  regularization  calculated  in  the  entire  image, 
(b)  Schematic  plot  showing  the  underestimation  of  the  displacement  (Data  +  reg.  curve)  with  unbiased  regularization  (refer  to  the  text),  (c),  (d)  The  calculated 
displacements  at  the  bottom  of  a  RF-line  at  2%  strain  and  6%  strain  levels  respectively  with  biased  and  unbiased  regularization  terms.  The  ground  truth  matches  the 
displacement  given  by  the  biased  regularization  almost  perfectly,  and  therefore  is  not  shown  in(c)  and  ( d)  not  to  block  the  biased  regularization  results.  The  length 
of  the  RF-line  is  2560  (49.3  mm),  (e)  The  SNR  values  corresponding  to  the  unbiased  regularization  calculated  by  omitting  the  bottom  300  samples  of  the  image, 
(f)  The  SNR  values  corresponding  to  the  biased  regularization  calculated  in  the  entire  image.  Note  that  the  scale  of  the  SNR  in  graph  (a)  is  much  smaller  than  that 
of  graphs  (e)  and  (f).  (a)  Unbiased  reg.  Entire  image,  (b)  Schematic  displacements,  (c)  Calculated  displacements  at  2%  strain,  (d)  Calculated  displacements  at  6% 
strain,  (e)  Unbiased  reg.  Top  of  the  image,  (f)  Biased  reg.  Entire  image. 


-  -  Based,  w  IRIS 

°0  2  4  6  8  1C 


stranS 

Fig.  4.  Lateral  strain  estimation  using  the  2D  AM  method  in  the  first  simulated 
phantom. 


blue  curve  increases  while  the  black  dashed  curve  decreases). 
Therefore,  one  cannot  conclude  from  this  experiment  that  higher 
a  is  beneficial  to  both  bias  and  variance.  To  prove  this,  we  de¬ 
signed  a  simulation  study  where  the  underlying  axial  strain  field 
continuously  varied  with  depth  and  the  lateral  and  out-of-plane 
strains  were  zero  (such  strain  field  is  not  physically  realizable). 
We  observed  that  the  absolute  value  of  the  bias  monotonically 
increases  with  a  with  both  unbiased  and  biased  regularizations. 
To  save  space,  we  do  not  present  the  full  results  here.  Similar 
curves  for  the  lateral  strain  field  is  shown  in  Fig.  6. 

The  second  simulation  experiment  is  designed  to  show  the 
effect  of  smoothness  weight  and  IRLS  threshold  CNR  when  the 
correlation  is  lower  in  parts  of  the  image  due  to  fluid  motion.  The 
phantom  contains  a  vein  oriented  perpendicular  to  the  image 
plane  (Fig.  7).  The  background  window  for  CNR  calculation  is 


(a)  (b) 


Fig.  5.  Bias  and  Variance  of  the  axial  strain  as  a  function  of  the  axial  regu¬ 
larization  weight  a.  The  ground  truth  axial  and  lateral  strain  fields  are  respec¬ 
tively  uniform  2%  and  2i/%  fields  ( v  —  0.49  is  the  Poisson’s  ratio).  The  solid 
blue  and  dashed  black  curves  both  correspond  to  unbiased  regularization  and 
the  solid  black  curve  corresponds  to  the  biased  regularization.  In  the  solid  blue 
and  solid  black  curves,  the  entire  image  is  included  in  the  calculation  of  the  bias 
and  noise.  In  the  dashed  black  curve  the  bottom  part  of  the  strain  field  which 
suffers  from  high  bias  [Fig.  3(b)]  is  excluded  from  the  calculation  of  the  bias 
and  noise.  ID  AM  and  2D  AM  have  very  similar  bias  and  variance.  The  curves 
with  and  without  IRLS  are  also  very  close.  Therefore  each  curve  corresponds  to 
ID  AM  or  2D  AM  with  or  without  IRLS.  (a)  Bias,  (b)  Variance. 


located  close  to  the  target  window  to  show  how  fast  the  strain  is 
allowed  to  vary,  a  property  related  to  the  spatial  resolution.  The 
maximum  CNR  with  IRLS  is  5.3  generated  at  T  —  0.005  and 
aa  =  38,  and  without  IRLS  is  4.8  at  aa  =  338.  Such  high  aa 
value  makes  the  share  of  the  data  term  in  the  cost  function  very 
small  and  causes  over- smoothing. 
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Fig.  6.  Bias  and  Variance  of  the  lateral  strain  as  a  function  of  the  axial  regularization  weight  a.  The  ground  truth  axial  and  lateral  strain  fields  are  respectively 
uniform  2%  and  2v%  fields  ( v  —  0.49  is  the  Poisson’s  ratio).  The  solid  blue  curve  corresponds  to  unbiased  regularization  and  the  dashed  and  solid  black  curves 
correspond  to  the  biased  regularization.  IRLS  is  not  used  in  the  solid  blue  and  dashed  black  curves,  (a)  Bias,  (b)  Variance. 


Fig.  7.  Measurements  in  (a)  are  in  mm.  In  (b),  a  scatterer  is  shown  in  the  bottom 
left  part  as  a  red  dot.  Its  displacement  is  calculated  by  interpolating  the  dis¬ 
placements  of  its  three  neighboring  nodes  on  the  mesh.  The  target  (circular) 
and  background  (rectangular)  windows  for  CNR  calculation  of  (d)  are  shown  in 

(c) .  (a)  Simulation  phantom,  (b)  Finite  element  mesh,  (c)  Finite  element  strain. 

(d)  CNR. 


A.  Displacement  Simulation 

To  study  the  performance  of  the  Kalman  filter,  we  simulate 
a  displacement  field  of  size  100  x  100  samples  whose  strain 
image  (calculated  using  least  squares  regression)  is  as  shown 
in  Fig.  8(a).  One  hundred  samples  in  the  axial  direction  corre¬ 
sponds  to  approximately  1.9  mm  (assuming  40  MHz  sampling 
rate),  and  100  samples  in  the  lateral  direction  corresponds  to 
10-25  mm  depending  on  the  probe.  To  be  consistent  with  the 
notations  of  Section  II-D,  let  c  j  denote  the  strain  values  of 
the  uncontaminated  image  in  (a).  We  then  contaminate  the  dis¬ 
placement  field  with  a  Gaussian  noise  with  standard  deviation 
of  1 .5  samples,  and  perform  least  squares  regression  to  calculate 
the  noisy  estimates  Zij  [Fig.  8(b)].  We  then  apply  the  Kalman 
filter  as  described  in  Section  II-D  to  the  noisy  estimates  in 
the  lateral  direction  (i.e.,  row-by-row).  The  posterior  estimates 
of  the  strain  values,  are  shown  in  (c).  The  strain  values  of 
the  shown  line  in  (a)-(c)  (at  i  =  50  samples)  is  shown  in  (d) 
and  (e)  [The  plot  in  (d)  around  the  step  in  magnified  in  (e)]. 
The  Kalman  filter  formulation  is  eliminating  the  noise  without 
over- smoothing  the  strain  image.  This  is  due  to  the  model  vari¬ 
ance  update  (27).  We  note  that  although  displacement  is  gen¬ 
erally  continuous  in  tissue,  its  spatial  derivation  (strain)  is  not: 
at  the  boundary  of  two  tissues  with  different  elasticity  moduli, 
strain  field  is  discontinuous. 


IV.  Experimental  Results 

For  experimental  evaluation,  RF  data  is  acquired  from  an 
Antares  Siemens  system  (Issaquah,  WA)  at  the  center  frequency 
of  6.67  MHz  with  a  VF10-5  linear  array  at  a  sampling  rate  of 
40  MHz.  Only  the  2D  AM  method  is  used  in  the  experimental 
results.  Phantom  results  and  patient  trials  are  presented  in  this 
section.  The  tunable  parameters  of  the  2D  AM  algorithm  are  set 
to  a  =  5,  j3a  =  10,  (3i  =  0.005  and  T  =  0.2  [(12)  and  (20)], 
and  the  tunable  parameters  of  the  DP  (run  for  the  seed  RF-line 
in  the  2D  AM  algorithm)  are  aa  —  ai  —  0.15  (1)  in  all  the 
phantom  results  (except  if  specified  otherwise).  In  the  patient 
results,  all  the  parameters  are  the  same  except  for  f3a  which  is 
increased  to  j3a  =  20  because  the  data  is  noisier.  The  strain  im¬ 
ages  in  all  the  patient  trials  are  obtained  using  the  least  squares 
regression  and  Kalman  filtering  as  described  in  Section  II-D. 

A.  Phantom  Results 

1)  Effect  of  Regularization  on  Residuals:  The  cost  function 
of  the  AM  method  (7)  is  composed  of  residuals  (i.e.,  the  data 
term)  and  the  regularization  terms.  The  AM  method  minimizes 
this  summation.  Therefore  the  AM  method  will  not  necessarily 
minimize  the  residuals.  We  now  show  that  the  data  term  alone 
is  nonconvex  and  has  many  local  minima.  Adding  the  regular¬ 
ization  term  will  eliminate  many  of  the  local  minima  and  makes 
optimization  of  the  data  term  easier.  This  is  in  addition  to  the  ef¬ 
fect  of  regularization  that  makes  the  displacement  field  smooth, 
a  generally  desired  attribute. 

The  effect  of  regularization  on  the  residuals  is  studied 
using  experimental  data.  An  elastography  phantom  (CIRS 
elastography  phantom,  Norfolk,  VA)  is  compressed  0.2  in 
axially  using  a  linear  stage,  resulting  in  an  average  strain  of 
6%.  Two  RF  frames  are  acquired  corresponding  to  before  and 
after  the  compression.  The  Young’s  elasticity  modulus  of  the 
background  and  the  lesion  under  compression  are  respectively 
33  kPa  and  56  kPa.  The  displacement  map  is  calculated  using 
the  2D  AM  method  and  the  residuals  corresponding  to  the 
displacement  map  are  obtained.  Fig.  9(a)-(c)  shows  the  axial 
and  lateral  strains  at  such  a  high  strain  rate  (minimum  of  2% 
and  maximum  of  11%).  The  mean  and  median  of  the  residuals 
in  the  entire  image  is  shown  in  (d).  One  could  expect  the 
graph  to  monotonically  increase  as  the  regularization  weight  a 
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Fig.  8.  (a)  shows  the  strain  field  calculated  using  least  squares  regression  of  the  uncontaminated  displacement  field,  (b)  depicts  the  strain  field  calculated  using  least 
squares  regression  of  the  contaminated  displacement  field,  (c)  shows  the  strain  field  calculated  from  the  noisy  measurements  of  (b)  using  the  proposed  Kalman  filter 
(KF  in  (b)  and  (c)  refers  to  Kalman  filter).  The  pixels  of  images  in  (a)  to  (c)  are  respectively  the  ground  truth  (unavailable)  strain  values  ,  the  noisy  measurements 
Zij,  and  posterior  strain  values  e<tj.  The  brightness  scale  in  (a)-(c)  is  the  same,  (d),  (e)  are  the  strain  estimation  at  the  horizontal  line  shown  in  (a)-(c).  (d)  is 
magnified  in  (e)  around  the  step.  The  Kalman  filter  removes  the  noise  while  keeping  the  image  sharp,  due  to  the  variable  model  noise  of  (27).  (a)  Ground  truth 
strain,  (b)  Strain  without  KF.  (c)  Strain  with  KF.  (d)  Strain  estimate,  (e)  Strain  estimates. 


increases,  since  the  difference  between  the  objective  function 
C  and  the  residuals  Y^ILx  p(ri)  *s  increased  as  a  is  increased. 
However,  the  residual  values  are  very  high  at  very  low  a.  There¬ 
fore,  numerical  minimization  of  YliLi  p(ri)  +  R(Ad)  gives 
a  smaller  value  for  Y^ILi  p(ri)  compared  to  trying  to  directly 
minimize  p(rd-  This  indicates  that  the  nonregularized 

cost  function  is  not  quasi-convex  and  is  very  hard  to  minimize. 

2)  Resolution  of  the  Strain  Images  Generated  With  AM:  The 
effect  of  the  regularization  on  spatial  resolution  is  evaluated  ex¬ 
perimentally  using  the  experimental  setup  of  the  previous  ex¬ 
periment.  The  compression  is  set  to  0.1  in  in  this  experiment. 
Fig.  10(a)  shows  the  strain  image  obtained  by  compression  the 
lesion  with  the  Young’s  modulus  of  56  kPa.  Spatial  resolution  is 
evaluated  using  modulation  transfer  function  (MTF),  an  estab¬ 
lished  method  for  estimating  the  spatial  resolution  of  medical 
imaging  systems  that  was  relatively  recently  extended  to  elas- 
tography  [58].  The  spatial  resolution  of  the  reconstructed  im¬ 
ages  is  determined  with  a  three-step  approach  [59],  [60].  First, 
the  edge  spread  function  is  computed  by  averaging  the  pixel 
values  across  the  background-inclusion  interface  [the  line  in 
Fig.  10(a)].  Second,  the  line  spread  function  (LSF)  is  computed 
by  differentiating  the  edge  spread  function.  Third,  the  MTF  is 
determined  by  computing  the  Fourier  transform  of  the  LSF  and 
normalizing  the  resulting  function  to  zero  spatial  frequency 

MTF  (k)  =  .  (29) 

Fig.  10(c)  shows  the  MTF  for  five  different  normalization  coeffi¬ 
cients  respectively.  Strain  results  are  obtained  with  a  regression 
window  of  length  2k  +  1  =  65  [Section  II-D].  Increasing  the 


regularization  weight  is  adversely  affecting  spatial  resolution. 
Spatial  resolution  is  defined  as  the  spatial  frequency  when  the 
value  of  MTF  is  0.1.  At  a  =  1,  a  —  2  and  a  —  4  this  value 
is  respectively  2  cycles/mm,  1  cycles/mm,  and  0.5  cycles/mm. 
In  addition  to  a ,  this  value  also  depends  on  the  length  of  the  re¬ 
gression  window  2  k  +  1. 

3)  Image  Quality  Versus  Axial  and  Lateral  Sampling  Rates 
of  the  RF-Data:  Sampling  rate  of  the  RF-data  usually  ranges 
from  20  to  50  MHz  depending  on  the  hardware  of  the  device. 
The  number  of  the  A-lines  provided  in  an  image  also  varies  sig¬ 
nificantly.  In  addition,  bandwidth  limitations  of  the  data  transfer 
can  impose  limits  on  the  size  of  the  image  for  real-time  oper¬ 
ations.  In  this  study,  we  downsample  the  RF-data  by  a  factor 
of  2-4  in  the  axial  direction  and  by  a  factor  of  2-8  in  the  lat¬ 
eral  direction.  Fig.  ll(a)-(g)  shows  axial  and  lateral  displace¬ 
ment  and  strain  images  of  the  CIRS  elastography  phantom  un¬ 
dergoing  maximum  axial  strain  of  5%.  Axial  sampling  rate  can 
be  reduced  by  a  factor  of  2  without  significant  impact  on  the 
strain  image  quality  [part  (h)].  Downsampling  the  images  in  the 
lateral  direction  by  a  factor  of  4  results  the  CNR  of  the  axial 
and  lateral  strain  images  to  drop  respectively  12%  (from  16.3 
to  14.3)  and  56%  (from  2.55  to  1.13)  as  shown  in  (i).  While  the 
axial  strain  is  robust  to  the  number  of  A-line  in  the  image  even 
at  a  high  strain  level  of  5%,  the  lateral  strain  is  sensitive  to  it 
(i).  Similar  study  with  lower  axial  strain  levels  shows  that  as  the 
axial  strain  decreases,  higher  downsampling  rates  in  both  axial 
and  lateral  directions  are  possible  without  a  large  impact  on  the 
results. 

4)  Kalman  Filter:  The  performance  of  the  Kalman  filter  is 
studied  using  the  RF-data  used  in  Fig.  9.  The  linear  least  squares 
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Fig.  9.  Phantom  experimental  results.  The  top  row  shows  axial  displacement  and  axial  strains  as  labeled  (KF  in  (c)  refers  to  Kalman  filter).  Average  axial  strain  and 
maximum  strain  are  approximately  6.6%  and  11%.  (d)  and  (e)  show  lateral  displacement  and  lateral  strain,  respectively,  (f)  shows  residuals  as  the  regularization 
weight  varies,  (a)  Axial  displacement  (mm),  (b)  Axial  strain,  (c)  Axial  strain  with  KF.  (d)  Lateral  displacement  (mm),  (e)  Lateral  strain,  (f)  Residuals. 


Fig.  10.  Phantom  experimental  results  showing  the  resolution  of  the  2D  AM.  (a)  Strain  image.  The  edge  spread  function  is  evaluated  along  the  vertical  line, 
(b)  The  strain  across  the  edge  [vertical  line  in  (a)]  for  the  five  shown  regularization  values,  (c)  The  MTF  calculated  across  the  vertical  line  in  (a).  Spatial  resolution 
is  defined  as  the  spatial  frequency  when  the  value  of  MTF  is  0. 1 .  (a)  Axial  strain,  (b)  Strain  profile,  (c)  MTF. 


differentiation  technique  is  applied  to  the  axial  displacement 
field  calculated  with  2D  AM,  resulting  in  [Fig.  12(a)].  The 
Kalman  filter  is  then  applied  to  Zij  measurements  of  (a),  giving 
the  posterior  measurements  of  (b).  Comparing  the  strain 
values  at  a  horizontal  line  of  (a)  and  (b),  the  noisy  z>ij  measure¬ 
ments  are  smoothed  in  the  lateral  direction  using  the  proposed 
Kalman  filter,  with  minimal  blurring  of  the  edge. 

B.  Clinical  Study 

Seven  patients  undergoing  open  surgical  radiofrequency  (RF) 
thermal  ablation  for  primary  or  secondary  liver  cancer  were  en¬ 
rolled  between  February  06,  2008  and  July  28,  2009.  All  pa¬ 
tients  enrolled  in  the  study  had  unresectable  disease  and  were 
candidates  for  RF  ablation  following  review  at  our  institutional 
multidisciplinary  conference.  Patients  with  cirrhosis  or  subop- 
timal  tumor  location  were  excluded  from  the  study.  All  patients 
provided  informed  consent  as  part  of  the  protocol,  which  was 


approved  by  the  institutional  review  board.  RF  ablation  was  ad¬ 
ministered  using  the  RITA  Model  1500  XRF  generator  (Rita 
Medical  Systems,  Fremont,  CA).  Strain  images  are  generated 
offline.  Some  preliminary  results  are  published  in  [15]. 

We  show  the  results  from  only  four  patients  due  to  space  lim¬ 
itations.  Fig.  13  shows  the  5-mode  scan,  the  strain  images  and 
CT  scans  performed  after  RF  ablation.  Tissue  is  simply  com¬ 
pressed  freehand  at  a  frequency  of  approximately  1  compres¬ 
sion  per  2  s  with  the  ultrasound  probe  without  any  attachment. 
The  shadow  in  Fig.  13(a)  at  20  mm  depth  is  produced  by  the 
thermal  lesion.  Note  that  it  is  not  possible  to  ascertain  the  size 
and  position  of  the  thermal  lesions  from  5-mode  images.  In  ad¬ 
dition,  the  thermal  lesion  has  different  appearances  in  the  three 
5-scans.  However,  the  thermal  lesions  show  very  well  as  hard 
lesions  in  the  strain  images.  After  gross  correlation  of  the  post 
ablation  CT  scan  and  the  thermal  lesion  in  the  strain  images,  the 
size  of  the  lesion  seems  to  correspond  well.  However,  a  more 
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Fig.  11.  Results  of  the  CIRS  elastography  phantom  at  5%  maximum  strain  at  different  axial  and  lateral  sampling  rates.  The  hard  lesion  is  spherical  and  has  a 
diameter  of  1  cm.  Downsampling  is  performed  by  simply  skipping  samples  in  the  axial  or  (and)  lateral  directions.  In  (c)  and  (f),  a  downsampling  ratio  of  2  is 
applied  in  both  axial  and  lateral  directions.  The  lateral  displacement  is  shown  in  number  of  samples  in  (d)-(f).  (h)  and  (i)  show  the  CNR  between  the  target  and 
background  windows  in  the  strain  images  as  the  axial  or  lateral  downsampling  rates  change.  The  target  and  background  windows  are  shown  in  the  axial  strain 
images  (a)-(c)  and  the  lateral  strain  image  (g).  In  (i),  the  lateral  strain  curve  is  not  calculated  for  downsampling  ratios  of  6  and  higher  because  the  background 
window  moves  out  of  the  image.  The  black  dashed  curve  with  the  highest  CNR  is  the  strain  obtained  with  the  Kalman  filter  (KF).  (a)  Axial  downsamp.  ratio  = 
2.  (b)  Lateral  downsamp.  ratio  =  2 .  (c)  Ax.-lat.  downsamp.  rat  =  2 .  (d)  Axial  downsamp.  ratio  =  2 .  (e)  Lateral  downsamp.  ratio  =  2 .  (f)  Ax.-lat.  downsamp. 
rat  =  2.  (g)  Lateral  downsamp.  ratio  =  2. 


rigorous  validation  of  the  size  and  shape  of  the  ablated  lesion  in 
the  elastography  image  is  underway  using  nonrigid  registration 
of  CT  and  ultrasound  images.  To  the  best  of  our  knowledge,  this 
is  also  the  first  demonstration  of  the  success  of  elastography  in 
imaging  the  thermal  lesion  in  an  in  vivo  human  experiment. 

We  have  also  acquired  patient  RF  data  of  liver  ablation  prior 
and  after  ablation  in  one  of  the  patient  trials.  Fig.  14  shows  the 
5-mode,  strain  and  venous  and  arterial  phase2  CT  images  ob¬ 
tained  before  ablation,  and  Fig.  15  shows  the  5-mode,  strain 
and  lateral  displacement  images  after  ablation.  In  Fig.  14,  the 
tumor  [marked  in  the  CT  images  (f)  and  (g)]  is  not  visible  in 

2CT  scans  are  performed  at  different  phases  after  intravenous  injection  of  a 
contrast  agent.  In  the  arterial  phase  (directly  after  injection  of  a  contrast  agent) 
arteries  will  enhance,  where  as  in  the  venous  phase  (30-60  s  after  injection)  the 
hepatic  parenchyma  and  veins  will  enhance. 


the  5-mode  image  (a),  but  is  clearly  visible  in  the  strain  images 
(b)  and  (c).  While  the  tissue  is  getting  compressed  with  the  ul¬ 
trasound  probe,  the  middle  hepatic  vein  (marked  as  5)  which 
is  only  4-8  cm  from  vena  cava  inferior  pulsates  at  high  ampli¬ 
tude.  The  graph  in  (e)  schematically  shows  the  probe  motion 
and  variations  in  the  diameter  of  the  vein.  Therefore,  the  vein 
can  look  soft  as  in  (c)  or  hard  as  in  (b)  depending  on  whether  its 
diameter  variation  is  in  the  same  [marked  by  ellipse  1  in  (e)]  or 
opposite  [marked  by  ellipse  2  in  (e)]  direction  as  the  probe  mo¬ 
tion.  The  effect  of  pulsation  of  vessels,  a  well-known  cause  of 
signal  decorrelation,  is  minimized  via  IRLS  resulting  in  a  low 
noise  strain  image.  In  addition,  since  the  2D  AM  method  gives  a 
dense  motion  field  (same  size  as  RF  data),  the  small  artery  at  the 
diameter  of  less  than  2  mm  [marked  as  4  in  (a)]  is  discernible  in 
(b)  from  the  low  pressure  portal  vein.  The  ablated  lesion  is  also 
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Fig.  12.  (a)  Shows  the  axial  strain  field  calculated  by  least  squares  regression  of  the  noisy  displacement  field,  (b)  depicts  the  strain  field  calculated  from  the  noisy 
measurements  of  (a)  using  the  proposed  Kalman  filter  (KF  in  (a)  and  (b)  refers  to  Kalman  filter).  The  pixels  of  images  in  (a)  and  (b)  are  respectively  the  least 
squares  measurements  zitJ ,  and  posterior  strain  values  e7  ^t .  (c)  shows  the  strain  estimation  at  the  17  mm  deep  horizontal  line  shown  in  (a)  to  (b).  The  Kalman  filter 
removes  the  noise  while  keeping  the  image  sharp,  due  to  the  variable  model  noise  of  (27).  (a)  Strain  without  KF.  (b)  Strain  with  KF.  (c)  Strain  plots. 
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Fig.  13.  In  vivo  images  of  the  thermal  lesion  produced  by  RF  ablation  therapy  of  liver  cancer.  All  images  acquired  after  ablation.  First,  second,  and  third  rows 
correspond  to  the  first,  second  and  third  patients  respectively.  The  thermal  lesion  shows  in  (b),  (f)  and  (j)  as  dark,  surrounded  by  normal  tissue  in  white.  The  lateral 
displacement  images  are  shown  in  number  of  samples  (they  do  not  immediately  carry  anatomical  information).  In  (b),  (d),  (f),  (h),  (j),  and  (1)  the  delineated  thermal 
lesions  is  shown.  The  nonunity  aspect  ratio  in  the  axes  of  the  5-mode  and  strain  images  should  be  considered  when  comparing  them  to  the  CT  scans,  (a)  5-mode 
patient  1.  (b)  Axial  strain,  (c)  Lateral  displacement,  (d)  CT  patient  1.  (e)  5-mode  patient  2.  (f)  Axial  strain,  (g)  Lateral  displacement,  (h)  CT  patient  2.  (i)  5-mode 
patient  3.  (j)  Axial  strain,  (k)  Lateral  displacement.  (1)  CT  patient  3. 


discernible  in  the  strain  images  of  Fig.  15(b)  and  (c).  We  believe 
the  soft  region  in  the  middle  of  the  two  hard  ablation  lesions  in 
(b)  and  (c)  (at  the  depth  of  25-30  mm  and  width  of  10-25  mm) 
is  not  close  to  any  of  the  10  tines  of  the  ablation  probe.  There¬ 
fore  because  of  its  proximity  to  veins  and  vessels  its  temperature 
has  remained  low. 


V.  Discussion 

The  resolution  of  the  method  is  formally  studied  in 
Section  IV-A  using  the  phantom  experiment.  Future  work 
will  include  more  intuitive  measures  for  resolution  in  terms 
of  the  smallest  detectable  target  as  a  function  of  its  elasticity 
difference  with  the  background. 
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Fig.  14.  In  vivo  images  of  the  fourth  patient  before  RF  ablation.  In  (a),  the  left  anterior  branch  of  portal  vein  is  marked  as  1  and  2  and  has  low  pressure  and 
therefore  compresses  easily.  Arteries  (marked  as  3  and  4)  and  the  middle  hepatic  vein  (marked  as  5)  however  pulsate  with  the  heart  beat  and  may  have  low  or  high 
pressure,  (b)  and  (c)  both  show  the  axial  strain  from  the  same  location  before  ablation.  They  are  calculated  at  two  different  phases  of  the  heart  beat.  The  cancer 
tumor  is  discernible  in  (b)  and  (c)  (regardless  of  the  systolic  or  diastolic  blood  pressure),  and  its  boundary  is  shown.  1  and  2  [as  marked  in  (a)]  correspond  to  the 
high  strain  area  in  both  (b)  and  (c).  Since  3,  4,  and  5  [as  marked  in  (a)]  pulsate,  they  may  look  hard  [as  in  (b)]  or  soft  [as  in  (c)].  (d)  shows  the  lateral  displacement. 
The  tumors  is  not  visible  in  this  image,  (e)  shows  the  motion  of  the  probe  and  the  variation  in  the  diameter  of  the  arteries  due  to  the  heart  beat  (refer  to  the  text), 
(f)  is  the  arterial  phase  and  (g)  is  the  venous  phase  contrast  CT  images.  The  numbers  1-5  mark  the  same  anatomy  as  (a),  (a)  5-mode  pre-ablation,  (b)  Axial  strain 
pre-ablation,  (c)  Axial  strain  pre-ablation,  (d)  Lateral  displacement  pre-ablation,  (e).  (f)  CT  pre-ablation,  (g)  CT  pre-ablation. 
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Fig.  15.  In  vivo  images  of  the  fourth  patient  after  RF  ablation.  Similar  to  Fig.  14,  the  hepatic  vein  (marked  as  5)  can  have  low  strain  [as  in  (b)]  or  high  strain  [as 
in  (c)]  values,  (a)  5-mode  post-ablation,  (b)  Axial  strain  post-ablation,  (c)  Axial  strain  postablation,  (d)  Lateral  displacement  postablation. 


The  cost  function  is  a  regularized  function  of  all  displace¬ 
ments  on  an  A-line.  This  makes  the  methods  robust  to  noise 
which  exist  throughout  the  image.  Besides,  the  AM  methods  are 
not  window-based  and  therefore  they  do  not  suffer  from  decor¬ 
relation  within  the  window.  As  a  result,  both  AM  methods  work 
for  strains  as  high  as  10%.  In  addition,  the  IRLS  outlier  rejec¬ 
tion  technique  makes  the  AM  methods  robust  to  local  sources 
of  decorrelation  such  as  out-of-plane  motion  of  movable  struc¬ 
tures  or  blood  flow. 

Global  stretching  assumes  a  constant  strain  across  the  depth 
and  stretches  one  of  the  RF-limes  accordingly.  It  is  shown  that  it 
enhances  the  quality  of  correlation  based  elastography  methods. 
The  reason  is  that  the  strain  of  each  point  can  be  assumed  to  be 
the  global  strain  (fixed  for  each  RF-line)  plus  some  perturbation, 
i.e.,  constant  strain  is  a  better  approximation  than  zero  strain. 
Biasing  the  regularization  is  motivated  by  the  same  reason  and 
involves  almost  no  additional  computational  cost. 


Improvement  in  the  SNR  and  CNR  achieved  with  Kalman 
filtering  differentiation  is  due  to  utilizing  the  (piecewise)  conti¬ 
nuity  of  the  strain  field.  One  could  think  of  a  unified  framework 
which  includes  both  the  2D  AM  and  the  Kalman  filtering  and 
directly  calculates  the  strain  field.  We  made  an  effort  to  formu¬ 
late  (15)  in  terms  of  strain  values.  Unfortunately,  the  coefficient 
matrix  in  the  left-hand  side,  became  a  full  matrix  for  our  de¬ 
sired  regularization.  Such  large  full  system  cannot  be  solved  in 
real-time. 

The  least  squares  differentiation  of  Section  II-D  can  be  in¬ 
corporated  in  the  Kalman  filter.  This  can  be  simply  done  by 
defining  the  state  at  each  point  to  be  the  displacement  and  the 
strain  of  that  point.  The  observed  variables  are  the  noisy  dis¬ 
placement  measurements  from  2D  AM.  Solving  for  the  state 
gives  a  strain  estimate  at  each  point.  However,  we  preferred 
to  follow  the  common  approach  of  first  finding  the  strain  by 
solving  least  squares.  In  addition,  the  axial  and  lateral  displace- 
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ments  can  be  considered  as  two  channels  of  a  measurement  and 
a  Kalman  filter  that  takes  into  account  both  intra-channel  (spa¬ 
tial)  and  inter-channel  variations  can  be  developed.  This  is  a  sub¬ 
ject  of  future  work. 

Lateral  displacement  estimation  with  2D  AM  is  of  order  of 
magnitude  less  accurate  than  the  axial  displacement  estimates. 
We  tested  the  following  algorithm  for  calculating  the  lateral  dis¬ 
placement  field  based  on  ID  AM:  run  ID  AM  to  find  the  axial 
displacement  field  A,  then  transpose  both  ultrasound  images  /i 
and  I2  and  run  ID  AM  again  using  A  calculated  in  the  previous 
step.  The  axial  displacement  field  calculated  for  the  transposed 
images  is  in  fact  the  lateral  displacement  of  the  original  im¬ 
ages.  Although  considerably  more  computationally  expensive 
than  2D  AM,  this  algorithm  did  not  improve  the  lateral  displace¬ 
ment  estimation.  Therefore  only  images  of  lateral  displacement 
are  provided  for  the  patient  trials  because  the  lateral  strain  did 
not  show  the  ablation  lesion.  This  is  in  accordance  with  recent 
work  [36]  which  only  shows  the  lateral  displacement.  A  2D  dis¬ 
placement  field  can  be  utilized  to  calculate  the  thermal  expan¬ 
sion  and  to  reconstruct  the  strain  tensor.  Incorporation  of  the 
synthetic  lateral  phase  [61]— [63],  into  2D  AM  to  further  improve 
the  accuracy  of  the  lateral  displacement  measurement  is  also  a 
subject  of  future  work. 

In  cases  where  the  two  ultrasound  frames  correlate  very 
poorly  throughout  the  image,  ID  AM  outperforms  2D  AM 
because  DP  is  run  for  the  entire  image  in  ID  AM.  However, 
in  those  cases  the  strain  images  are  of  very  low  quality  even 
with  ID  AM.  In  cases  where  the  images  correlate  reasonably, 
the  2D  AM  algorithm  slightly  outperforms  ID  AM  in  terms 
of  the  SNR  of  the  axial  strain  as  shown  in  Fig.  3(e)  and  (f). 
Also,  ID  AM  and  2D  AM  are  very  similar  in  terms  of  bias  and 
variance  as  mentioned  in  the  caption  of  the  Fig.  5.  And  finally, 
2D  AM  is  more  than  10  times  faster  than  ID  AM  because  it 
eliminates  the  redundant  calculations  in  the  DP  step  of  ID  AM. 
This  is  important  considering  that  there  are  combinatorial  many 
ways  of  choosing  two  frames  for  elastography  from  a  sequence 
of  images.  Having  a  fast  algorithm,  like  2D  AM,  makes  it 
plausible  to  invest  time  to  perform  real-time  frame  selection, 
an  area  that  we  are  currently  working  on  [16],  [64]. 

Recent  work  [65]  has  attempted  to  reconstruct  elasticity  from 
the  displacement  field  for  monitoring  thermal  ablation.  It  has 
also  shown  that  [66]  compared  to  strain  images,  elasticity  im¬ 
ages  have  both  higher  correlation  with  the  ablation  zone  and 
give  higher  CNR.  Another  work  [67]  has  utilized  the  solution 
of  the  elasticity  reconstruction  to  improve  motion  estimation  in 
an  iterative  framework.  Calculation  of  the  elasticity  modulus  in 
our  ablation  monitoring  trials  is  an  area  of  future  work. 

Statistical  analysis  of  the  residuals  is  a  subject  of  future  work. 
The  sum  of  squared  differences  used  as  the  similarity  metric  in 
our  cost  function  is  suitable  if  ultrasound  noise  can  be  modeled 
as  additive  Gaussian  noise.  However,  ultrasound  noise  is  not 
simply  additive  Gaussian  and  it  has  been  shown  that  similarity 
metrics  that  model  the  noise  process  considering  physics  of  ul¬ 
trasound  give  more  accurate  results  [68].  Performance  of  the  2D 
AM  method  for  images  that  are  not  fully  developed  speckles 
(i.e.,  have  few  scatterers  per  resolution  cell)  is  also  a  subject  of 
future  work. 

Current  implementations  of  the  ID  AM  and  2D  AM  take, 
respectively,  0.4  s  and  0.04  s  to  generate  strain  images  (axial 
for  ID  AM  and  axial  and  lateral  for  2D  AM)  of  size  1000  x  100 


on  a  3.8  GHz  P4  CPU.  DP  contributes  to  more  than  90%  of  the 
running  time  of  the  ID  AM,  and  that’s  why  it  is  slower  than  2D 
AM  where  DP  is  only  run  for  a  single  A-line.  The  running  time 
of  both  methods  changes  linearly  with  the  size  of  the  image. 

VI.  Conclusion 

Two  regularized  elastography  methods,  ID  AM  and  2D  AM, 
are  introduced  for  calculating  the  motion  field  between  two  ul¬ 
trasound  images.  They  both  give  dense  subsample  motion  fields 
(ID  AM  gives  subsample  axial  and  integer  sample  lateral  and 
2D  AM  gives  subsample  axial  and  lateral)  in  real-time.  The  size 
of  the  motion  fields  is  the  same  as  the  size  of  the  RF-data  (ex¬ 
cept  for  few  samples  from  the  boundary  whose  displacements 
are  not  calculated).  Such  dense  motion  fields  lead  to  dense  strain 
fields  which  are  critical  in  detecting  small  lesions.  The  prior  of 
tissue  motion  continuity  is  exploited  in  the  AM  methods  to  min¬ 
imize  the  effect  of  signal  decorrelation.  The  regularization  term 
is  biased  with  the  average  strain  in  the  image  to  minimize  un¬ 
derestimation  of  the  strain  values.  Parts  of  the  image  that  have 
very  low  correlation  are  treated  as  outliers  and  their  effect  is 
minimized  via  IRLS.  The  strain  image  is  calculated  by  differ¬ 
entiating  the  motion  fields  using  least  squares  regression  and 
Kalman  filtering.  The  performance  of  the  proposed  elastography 
algorithms  is  analyzed  using  Field  II  and  finite  element  simula¬ 
tions,  and  phantom  experiments.  Clinical  trials  of  monitoring 
RF  ablation  therapy  for  liver  cancer  in  four  patients  are  also 
presented.  An  implementation  of  the  2D  AM  method,  the  least 
squares  regression  and  the  Kalman  filter  in  MATLAB  mex  func¬ 
tions,  as  well  as  some  of  the  phantom  and  patient  RF  data  used  in 
this  work  are  available  for  academic  research  and  can  be  down¬ 
loaded.3 
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Abstract.  The  clinical  feasibility  of  2D  elastography  methods  is  hin¬ 
dered  by  the  requirement  that  the  operator  avoid  out-of-plane  motion  of 
the  ultrasound  image  during  palpation,  and  also  by  the  lack  of  volumet¬ 
ric  elastography  measurements.  In  this  paper,  we  develop  and  evaluate 
a  3D  elastography  method  operating  on  volumetric  data  acquired  from 
a  3D  probe.  Our  method  is  based  on  minimizing  a  cost  function  using 
dynamic  programming  (DP).  The  cost  function  incorporates  similarity 
of  echo  amplitudes  and  displacement  continuity.  We  present,  to  the  best 
of  our  knowledge,  the  first  in-vivo  patient  studies  of  monitoring  liver 
ablation  with  freehand  DP  elastography.  The  thermal  lesion  was  not  dis- 
cernable  in  the  B-mode  image  but  it  was  clearly  visible  in  the  strain 
image  as  well  as  in  validation  CT.  We  also  present  3D  strain  images 
from  thermal  lesions  in  ex-vivo  ablation.  Good  agreement  was  observed 
between  strain  images,  CT  and  gross  pathology. 


1  Introduction 

Hepatocellular  carcinoma  (HCC)  is  one  of  the  most  common  tumors,  caus¬ 
ing  662,000  deaths  worldwide  annually.  Minimally  invasive  RF  ablation  [1]  has 
gained  much  interest  recently  since  only  10%  to  20%  of  patients  with  HCC  are 
amenable  to  traditional  therapy  of  surgical  resection  of  the  tumor.  In  RF  ab¬ 
lation,  an  electrode  is  placed  into  the  tumor  to  cauterize  it  [1].  Monitoring  the 
ablation  process  in  order  to  document  adequacy  of  margins  during  treatment 
is  a  significant  importance.  Ultrasonography  is  the  most  common  modality  for 
both  target  imaging  and  for  ablation  monitoring.  However,  ultrasonographic  ap¬ 
pearance  of  ablated  tumors  only  reveals  hyperechoic  areas  due  to  microbubbles 
and  outgasing  and  cannot  adequately  visualize  the  margin  of  tissue  coagulation. 

Accordingly,  ultrasound  elastography  (Ophir  et  al,  1991)  has  emerged  as  a 
useful  augmentation  to  conventional  ultrasound  imaging.  Elastography  has  been 
used  for  monitoring  RF  ablation  [2],  [3]  by  observing  that  ablated  region  is 
harder  than  surrounding  tissue.  In  the  most  common  variation  of  elastography, 
ultrasound  images  are  captured  while  the  tissue  is  being  compressed,  and  images 
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are  processed  to  provide  a  grid  of  local  displacement  measurements.  These  dis¬ 
placement  fields  are  then  used  to  determine  the  elastic  properties  of  the  tissue 
at  each  grid  location.  The  grid  of  calculated  elastic  properties  can  be  displayed 
as  an  image. 

Elastography  is  computationally  expensive,  making  it  challenging  to  display 
strain  images  in  real-time.  Real-time  feedback,  however,  is  required  for  image 
guided  ablation  operations.  Another  aspect  is  that  signal  decorrelation  between 
the  pre-  and  post-compression  images  induces  significant  noise  in  the  obtained 
displacement  map  and  is  one  of  the  major  limiting  factors  in  elastography  [4]. 
Methods  based  on  cross-correlation  and  phase  zero  estimation  are  currently  the 
most  popular  real-time  elastography  techniques  which  provide  fast  and  accurate 
motion  tracking.  In  RF  ablation,  however,  high  decorrelation  between  pre-  and 
post-compression  images  results  in  high  noise  in  the  strain  images  obtained  using 
cross-correlation  [3].  Phase  zero  estimation  methods  require  an  estimate  of  the 
center  frequency  of  the  ultrasound  RF  signal,  which  varies  with  depth  due  to 
frequency-dependent  attenuation  in  tissue  [5].  This  variation  can  be  significant 
in  RF  ablation,  leading  to  poor  displacement  estimation  [5] . 

We  have  recently  developed  a  real-time  2D  elastography  method  based  on 
dynamic  programming  (DP)  [6].  The  method  is  more  robust  to  signal  decorrela¬ 
tion  than  standard  cross-correlation  methods  and  is  therefore  a  good  candidate 
for  ablation  monitoring  where  being  real-time  and  robustness  to  noise  are  criti¬ 
cally  important.  Here  we  report,  to  the  best  of  our  knowledge,  the  first  in-vivo 
patient  results  on  monitoring  RF  ablation  with  2D  DP  elastography  and  corrob¬ 
orate  the  results  with  CT  scans.  As  initial  clinical  studies  revealed  limitations 
of  2D  elastography  in  monitoring  the  thermal  ablation,  we  were  compelled  to 
progress  toward  3D.  We  think  the  readership  will  find  it  informative  to  see  how 
our  concept  and  methodology  evolved.  We  extend  our  DP  method  to  operate  on 
3D  volumes.  The  benefits  of  3D  strain  imaging  of  RF  ablation  are  two- fold:  1)  3D 
imaging  eliminates  the  need  to  image  the  same  plane  while  palpating  the  tissue, 
which  can  be  very  difficult  in  the  presence  of  breathing  and  cardiac  motion,  and 
2)  3D  imaging  allows  more  precise  monitoring  of  temperature  deposition  which 
exhibits  variations  in  3D,  particularly  in  the  presence  of  blood  vessels  which  act 
as  heat  sinks.  Previous  work  has  generated  3D  elastography  by  moving  a  con¬ 
ventional  2D  probe  out-of-plane  using  mechanical  guidance  [7,8]  or  freehand  [9]. 
In  recent  work  by  Treece  et  al.  [10]  and  Fisher  et  al.  [11]  a  3D  probe  is  used 
to  acquire  3D  elastogrophy,  using  phase  zero  and  cross-correlation  based  motion 
tracking  methods  respectively.  Here,  we  use  3D  probe  to  acquire  3D  data  and 
introduce  a  3D  DP  motion  tracking  algorithm.  We  show  that  3D  elastography 
can  be  successfully  used  to  monitor  ablation  in  3D. 

2  3D  Displacement  Estimation  Using  DP 

Compared  to  other  optimization  techniques,  DP  is  an  efficient  non-iterative 
method  of  global  optimization  [12,13].  We  have  recently  developed  a  real-time 


460 


H.  Rivaz  et  al. 


2D  elastography  method  using  DP  [6] .  In  DP  elastography,  a  cost  function  which 
incorporates  similarity  of  echo  amplitudes  and  displacement  continuity  is  mini¬ 
mized.  Since  data  alone  can  be  insufficient  to  solve  ambiguities  of  motion  track¬ 
ing  due  to  signal  decorrelation,  the  physical  priors  of  tissue  motion  continuity 
increases  the  robustness  of  the  technique  [6].  We  have  showed  that  DP  gener¬ 
ates  high  quality  strain  images  of  freehand  palpation  elastography  with  up  to 
10%  compression,  indicating  that  the  method  is  more  robust  to  signal  decor¬ 
relation  (caused  by  scatterer  motion  in  high  axial  compression  and  non-axial 
motions  of  the  probe)  in  comparison  to  the  standard  correlation  techniques. 
The  method  operates  in  less  than  lsec  and  is  thus  also  suitable  for  real  time 
elastography. 

Here,  we  extend  this  method  to  operate  on  3D  volumes.  Devising  a  DP  algo¬ 
rithm  for  optimization  involves: 

1.  Breaking  the  total  optimization  cost  into  a  sum  of  individual  costs,  such  that 
each  cost  corresponds  to  a  discrete  decision.  The  decisions  should  follow  each 
other  sequentially  and  the  cost  corresponding  to  each  decision  should  only 
depend  on  the  previous  and  not  the  future  decisions  (causality). 

2.  Determining  what  decisions  are  possible  at  each  stage. 

3.  Writing  a  recursion  on  the  optimal  cost  from  the  first  stage  to  the  final  stage. 

Let  gj(i)  be  the  intensity  of  the  ith  sample  (axial  direction),  jth  A-line  (lateral 
direction)  and  kth  frame  (out-of-plane  direction)  of  the  pre-compression  ultra¬ 
sound  volume.  Let  g^di  +  da)  correspond  to  the  post-compression  volume 
where  da,  di  and  de  represents  axial,  lateral  and  elevational  displacements  re¬ 
spectively,  and  the  size  of  the  volume  be  m  x  n  x  p.  The  difference  between  the 
two  signals,  A,  can  be  quantified  using  sum  of  absolute  differences  (SAD),  which 
is  computationally  inexpensive  to  compute  and  has  been  shown  to  have  good 
robustness  against  outliers  [14]: 


A(i,j,  k,  da,  dh  de)  =  g%(i)  -  (i  +  da) 


(1) 


where  the  axial,  lateral  and  elevational  search  ranges  are  limited  by  da?m^n  < 


is  the  smoothness  regularization.  The  cost  function  at  each  point  i,  j  and  k  is 


Cj(da,di,de,i )  =  A(i,j,k,da,di,de )  +  wiR(da,  dt,  de,  d 
C^a,  Sh6e,  i-l)  +  C^-liSg,  5u  6e,  i) 


+  w2R(da,  di,  de,  6a,  Si,5e) 


+  min 

8a,8i,5( 


2 
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where  w\  is  a  weight  for  governing  smoothness  in  the  elevational  direction 
and  W2  is  a  weight  for  governing  axial  and  lateral  smoothness1.  Generally,  the 
optimum  values  of  Sa,Si,Se  should  be  sought  in  the  entire  da:Tnax\  x 

[di,min  di^max]  x  [de,min  d^max\  space.  However,  since  the  strain  value  is  low 
in  elastography,  it  is  expected  and  desired  that  at  each  sample  of  RF  data, 
the  change  between  the  displacement  of  a  sample  and  its  previous  sample  is 
not  more  than  1.  Therefore,  the  search  range  is  limited  to  the  nine  values  of 
{da  -  1,  da,  da  +  1}  x{di  —  1,  di ,  di  +  1}  x{de  —  1,  de,  de  +  1},  which  results 
in  a  significant  gain  in  speed.  This  limit  on  the  search  range  does  not  affect  the 
results  even  in  a  high  strain  of  10%:  Ad  is  zero  for  nine  samples  and  one  for 
the  tenth  sample  on  average.  For  memoization  [12],  £a,  Si  and  Se  values  that 
minimize  the  cost  function  are  stored: 


Mj(i,  di, di,de)  = 


(4) 


arg  mm  . 

Sa,S„Se  1 


(Ck(6a,5l,5e,i~  1)  +  C*  i(*a,  Sh5e,  i) 

m  <  — - - - 


W2R{da-,  di ,  de,  Sa ,  Si ,  8e) 


The  cost  function  Cj  is  calculated  for  i  =  1  •  •  •  m,  da  =  da,mm  *  *  *  ^a,macc5  di  = 
di,min  •  •  •  di,max  and  de  =  de^min  •  •  •  de> max.  The  minimum  cost  at  i  =  m  gives  the 
displacement  of  this  point,  which  is  traced  back  to  i  =  1  using  the  M  function  to 
calculate  the  three  axial,  lateral  and  elevational  displacements  ( D  =  (da,  d/,  de))\ 

DkAi)=  arg  min  {Ck{da,  dt,  de,  i)}  ,  i  =  m 

J  da,di,de  J 

D$(i)  =  M(i  +  l,D$(i  +  l)),  i  =  1  ■  ■  -m  —  1  (5) 

This  gives  all  three  displacements  simultaneously,  in  contrast  with  other  3D  elas¬ 
tography  methods  which  give  displacement  in  each  direction  in  separate  steps. 

Further  speed-up  is  achieved  by  downsampling  the  signal  g(i)  in  the  axial 
direction  by  a  factor  of  f3  to  and  comparing  it  with  the  unaltered  signal 

g'(i).  This  is  done  by  simply  skipping  [3—1  samples  from  g(i)  and  performing  DP 
on  the  [3th  sample  as  illustrated  in  Figure  1  left.  This  generates  integer  displace¬ 
ment  estimations  at  m/ (3  samples.  The  displacement  of  the  skipped  samples  is 
then  simply  approximated  by  the  linear  interpolation  of  two  neighboring  points 
whose  displacements  are  calculated,  as  an  initial  guess  for  the  next  step. 

The  displacement  estimates  are  then  refined  to  subpixel  displacement  esti¬ 
mation  in  the  axial  direction.  The  original  signal  g(i)  (not  downsampled)  is 
compared  with  g'(i  +  d)  upsampled  by  a  factor  of  7  (Figure  1  right)  in  the  axial 

1  The  inclusion  of  the  cost  of  the  previous  line  (Cj_1(-  •  •))  guarantees  lateral  smooth¬ 
ness.  Instead,  we  could  force  the  displacements  of  each  pixel  to  be  similar  to  the 
displacements  of  the  neighboring  pixel  in  the  previous  A-line,  similar  to  what  we  did 
in  the  wiR(-  •  •)  term  to  enforce  elevational  smoothness.  The  former  is  preferred  since 
a  wrong  displacement  estimation  does  not  affect  the  neighboring  A-line’s  displace¬ 
ment  estimation.  However,  it  requires  the  Cf^_1(*  •  •)  to  be  kept  until  the  calculation 
of  Cj( •  •  •)  is  completed.  Therefore,  at  each  time  only  two  cost  finctions  are  stored  in 
the  memory,  making  the  memory  requirement  independent  of  the  number  of  A-lines. 


462 


H.  Rivaz  et  al. 


Fig.  1.  In  the  left,  the  cost  function  C  is  shown  when  DP  is  performed  on  g*(i)  (g(i) 
downsampled  by  a  factor  of  /3)  and  g'(i)  (not  downsampled).  Hashed  squares  indicate 
no  cost  calculation  is  performed  due  to  downsampling  of  g(i),  and  white  and  black 
representing  low  and  high  cost  values  respectively.  Displacement  is  calculated  at  m//3 
samples  in  this  stage  (/3  =  3  in  this  figure).  In  right,  a  new  cost  function  around  the 
optimum  path  of  the  first  stage’s  cost  function  (the  dashed  line)  is  created,  giving  a 
l/y  =?  1/2  pixel  displacement  accuracy  at  m  samples. 


direction  using  parabolic  interpolation.  Repeating  the  refinement  procedure  n 
times  results  in  a  refinement  factor  of  l/y72.  The  code  runs  in  approximately 
30sec  for  a  typical  volume  on  a  3.8GHz  P4  CPU. 

In  cross  correlation  methods,  subsample  displacement  is  usually  achieved  by 
interpolation  of  the  correlation  function,  which  is  subject  to  bias  and  jitter  [15]. 
Here  we  interpolate  the  original  RF  data  instead,  which  is  shown  to  have  similar 
performance  [15].  Although  cosine-fit  outperforms  parabolic-fit  interpolation  in 
terms  of  bias  and  jitter  [15],  the  latter  is  used  here  for  computational  simplicity. 

3  Results 

We  first  present  in-vivo  elastographic  monitoring  of  RF  ablation  therapy  of  HCC 
in  human  during  surgery  using  2D  DP  elastography.  RF  ablation  was  adminis¬ 
tered  using  the  RITA  Model  1500X  RF  generator  (Rita  Medical  Systems,  Fre¬ 
mont,  CA).  Ultrasound  RF  data  is  acquired  from  an  Antares  Siemens  system 
(Issaquah,  WA)  with  a  7.27MHz  linear  array  at  a  sampling  rate  of  40MHz.  Fig¬ 
ure  2  shows  the  B-mode  scan,  the  strain  image  obtained  using  the  DP  method 
and  CT  scans  performed  after  RF  ablation  (first  and  second  row  corresponding 
to  first  and  second  patient  respectively).  Tissue  is  simply  compressed  freehand 
with  the  ultrasound  probe  without  any  attachment.  The  shadow  in  Figure  2(a) 
at  20mm  depth  is  produced  by  the  thermal  lesion.  Note  that  it  is  not  possible 
to  ascertain  the  size  and  position  of  the  thermal  lesions  from  B-mode  images.  In 
addition,  the  thermal  lesion  has  different  appearances  in  the  two  B-scans.  How¬ 
ever,  the  thermal  lesions  show  very  well  as  hard  lesions  in  both  the  strain  images. 
The  size  of  the  thermal  lesion  in  the  strain  images  and  in  CT  scans  are  also  in 
accordance.  The  strain  images  provide  with  higher  contrast  of  the  thermal  le¬ 
sion  and  lower  noise  in  the  image,  compared  to  the  strain  images  of  RF  ablation 
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Fig.  2.  In-vivo  images  of  the  thermal  lesion  produced  by  RF  ablation  therapy  of  HCC, 
first  and  second  row  corresponding  to  the  first  and  second  patients,  (a)  V  (d)  13- 
scan  after  RF  ablation.  The  shadow  in  (a)  indicates  the  presence  of  thermal  lesion. 
It  is  almost  impossible  to  ascertain  the  size  and  position  of  the  thermal  lesion  from 
the  B-scans.  (b)  &;  (e)  Strain  images  after  RFA  ablation,  generated  using  2D  DP 
elastography  and  freehand  palpation  of  the  liver  tissue.  The  thermal  lesion  is  visible  in 
dark  surrounded  by  normal  tissue  in  white,  (c)  &;  (f)  Post-ablation  CT  scans,  with  the 
delineated  thermal  lesions  (The  non-unity  aspect  ratio  in  the  axes  of  the  B-mode  and 
strain  images  should  be  considered  when  comparing  them  with  the  CT  scans). 


reported  in  [3,8]  which  are  obtained  with  cross-correlation.  To  the  best  of  our 
knowledge,  this  is  also  the  first  demonstration  of  the  success  of  elastography  in 
imaging  the  thermal  lesion  in  an  in-vivo  human  experiment. 

A  Radionics  device  (Valleylab,  Boulder,  CO)  is  used  for  ex-vivo  RF  ablation. 
For  3D  elastography,  we  use  a  3D  probe  that  consists  of  a  curvilinear  array  that 
is  mechanically  rotated  to  scan  a  volume.  Ultrasound  RF  data  is  acquired  from 
an  Ultrasonix  system  (Vancouver,  BC)  at  4.5MHz  frequency,  20MHz  sampling 
rate  and  30%  bandwidth.  Figure  3  shows  the  experimental  setup  and  results. 
Comparing  strain  images  obtained  during  ablation  and  after  ablation,  the  growth 
of  the  thermal  lesion  can  be  observed.  There  is  also  a  good  agreement  between 
the  size  of  the  lesion  in  axial  and  lateral  directions  in  the  post-ablation  strain 
images  and  gross  pathology  photograph.  The  ablation  goes  beyond  the  field  of 
view  of  the  3D  ultrasound  probe  in  the  elevational  direction.  The  volumetric 
elastography  contrast  to  noise  ratio  (CNR)  [2]  (CNR  =  where  st  and 

Sb  are  the  spatial  strain  average  of  the  target  and  background  and  of  and  of  are 
the  spatial  strain  variance  of  the  target  and  background)  between  two  3mm  x 
3mm  x  3mm  cubes,  one  in  the  thermal  lesion  and  the  other  half  way  between  the 
liver  surface  and  the  lesion,  is  3.4.  Note  that  due  to  the  lateral  and  elevational 
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Fig.  3.  Ex- vivo  liver  RF  ablation  experiment.  The  ablation  power  is  set  to  8W  for 
10  min  and  the  cooler  is  turned  off  throughout  the  experiment,  (a)  The  experimental 
setup.  The  passive  arm  is  holding  the  3D  probe  and  the  liver  is  contained  in  the  gelatin 
phantom,  (b)  The  liver  sample  after  ablation  (cut  into  four  pieces)  with  the  thermal 
lesion,  (c)-(f)  3D  strain  images  6min  after  the  start  of  the  ablation  (target  temperature 
reached  90° C  in  3min  and  was  approximately  constant  in  the  next  7min  of  ablation). 
(g)-(j)  3D  strain  images  after  the  ablation,  (k)-(n)  3D  B-mode  images  after  the  ablation. 
Each  650  axial  samples  correspond  to  25mm  in  the  strain  and  B-mode  images.  Each 
40  lateral  samples  correspond  to  16.6mm  and  30.4mm  on  the  top  and  bottom  of  each 
image  (different  values  on  the  top  and  bottom  since  the  probe  is  curved).  Each  80 
elevational  samples  correspond  to  3.4mm  and  6.3mm  on  the  top  and  bottom  of  each 
image. 


regularization,  DP  elastography  is  working  reasonably  in  the  presence  of  the 
ablation  needle. 
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4  Discussion  and  Conclusion 

Previous  work  has  shown  promise  in  monitoring  ablation  in  2D  ex-vivo  and  ani¬ 
mal  experiments.  In  this  paper,  we  present  high  quality  in-vivo  2D  strain  images 
of  thermal  lesions  and  compared  them  to  post-ablation  CT  data.  Comparison  is 
more  qualitative,  however,  since  strain  images  are  2D  and  CT  data  is  3D  and 
ultrasound  is  not  tracked.  We  also  present  formulation  and  experimental  results 
of  a  3D  strain  imaging  system  based  on  DP.  In  DP,  we  regularize  the  problem  of 
3D  displacement  estimation:  regularization  in  2D  is  shown  to  increase  robustness 
[6].  As  a  result,  no  post  processing  step  such  as  median  filtering  is  performed. 

We  demonstrate  the  feasibility  of  3D  elastography  monitoring  of  RF  ablation 
for  the  first  time  using  a  3D  probe;  however,  we  are  planning  for  a  comprehensive 
comparison  of  the  3D  DP  with  other  3D  strain  imaging  techniques  [10,11].  The 
lateral  and  elevational  search  is  performed  only  to  increase  the  quality  of  the 
axial  strain:  the  lateral  and  elevational  displacements  are  integer  values  and  are 
not  suitable  for  calculating  strain.  Good  volumetric  CNR  between  the  thermal 
lesion  and  background  suggests  that  the  regularization  is  not  adversely  affecting 
CNR.  However,  a  study  similar  to  [6]  on  the  effect  of  the  3D  regularization  on 
the  CNR  and  resolution  should  be  done.  Having  an  elastography  system  for  3D 
ablation  monitoring  with  promising  ex-vivo  results,  in-vivo  patient  studies  under 
our  active  Institutional  Review  Board  (IRB)  approval  are  to  commence. 
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ABSTRACT 

Out-of-plane  motion  in  freehand  3D  ultrasound  can  be  estimated  using  the  correlation  of  corresponding  patches, 
leading  to  sensor  less  freehand  3D  ultrasound  systems.  The  correlation  between  two  images  is  related  to  their 
distance  by  calibrating  the  ultrasound  probe:  the  probe  is  moved  with  an  accurate  stage  (or  with  a  robot  in 
this  work)  and  images  of  a  phantom  are  collected,  such  that  the  position  of  each  image  is  known.  Since  parts 
of  the  calibration  curve  with  higher  derivative  gives  lower  displacement  estimation  error,  previous  work  limits 
displacement  estimation  to  parts  with  maximum  derivative.  In  this  paper,  we  first  propose  a  novel  method  for 
exploiting  the  entire  calibration  curve  by  using  a  maximum  likelihood  estimator  (MLE).  We  then  propose  for 
the  first  time  using  constrains  inside  the  image  to  enhance  the  accuracy  of  out-of-plane  motion  estimation.  We 
specifically  use  continuity  constraint  of  a  needle  to  reduce  the  variance  of  the  estimated  out-of-plane  motion. 
Simulation  and  real  tissue  experimental  results  are  presented. 

Keywords:  3D  ultrasound,  Speckle  decorrelation,  Fully  developed  speckle 


1.  INTRODUCTION 

Most  common  techniques  for  acquiring  3D  ultrasound  data  are  oscillating  head  probes  and  freehand  3D  ultra¬ 
sound.  In  oscillating  head  probes,  a  ID  ultrasound  transducer  is  automatically  swept  inside  the  probe,  enabling 
3D  image  acquisition.  In  freehand  3D  ultrasound,  a  position  sensor  is  attached  to  an  ordinary  probe  which  is 
swept  over  the  desire  region  by  the  clinician. 

Freehand  3D  ultrasound  is  inexpensive,  works  with  the  existing  2D  probes,  and  allows  arbitrary  3D  volume 
acquisition.  However,  the  need  for  the  additional  sensor  makes  it  difficult  to  use.  Sensorless  volume  reconstruction 
of  freehand  3D  ultrasound  is  possible  using  the  information  in  the  images  themselves:  out  of  plane  motion 
estimation  can  be  obtained  from  image  correlation,1  which  is  the  focus  of  this  work,  while  in  plane  motion  can 
be  estimated  through  image  registration2-4  or  by  using  techniques  similar  to  elastography.5-7 

The  granular  appearance  of  ultrasound  images  is  the  key  factor  in  out-of-plane  motion  estimation  (Figure  1). 
Each  pixel  in  an  ultrasound  image  is  formed  by  the  back-scattered  echoes  from  an  approximately  ellipsoidal  region 
called  the  resolution  cell.8  The  interference  of  scatterers  in  a  resolution  cell  creates  the  granular  appearance  of 
the  ultrasound  image,  called  speckle.  Although  of  random  appearance,  speckle  pattern  is  identical  if  the  same 
object  is  scanned  from  the  same  direction  and  under  the  same  focusing  and  frequency.  Speckle  characterization 
is  essential  in  many  areas  of  quantitative  ultrasound.  In  this  work,  it  is  a  prerequisite  for  speckle-based  distance 
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estimation.  We  use  low  order  moments  to  discriminate  fully  developed  speckle  (FDS)  patches  versus  coherent 
speckle  patches.9 


R  =  SNR  = 


(^) 

\l{A2vr)  -  (^)2 


(1) 


S  =  skewness  = 


((Av°  -  ( Av *))3) 
(( A2Vs )  —  ( AVs)2)i 


(2) 


where  A  is  the  amplitude  of  the  ultrasound  RF  envelope  in  the  analysis  patch,  vr  and  vs  are  the  signal  powers 
and  (•  •  •)  denotes  the  mean.  Here  we  use10  vr  =  2vs  =  1.  An  elliptical  discrimination  function  is  calculated  in 
the  R-S  plane  by  performing  principal  component  analysis  (PCA)  on  the  data  from  simulated  FDS  patches.10 
A  patch  is  then  classified  as  FDS  if  its  R-S  duple  falls  inside  this  ellipse. 


Having  found  FDS  patches  in  two  ultrasound  images,  the  correlation  between  them  is  used  for  estimating  the 
distance  between  the  two  images.11  The  R-S  metric  requires  approximately  3500  pixels  per  patch  (depending 
on  the  correlation  of  data12),  but  such  large  patches  (which  are  rectangles)  of  FDS  are  unlikely  to  be  found  in 
real  tissue  because  of  its  inhomogeneity.11  Gee  et  al.13  proposed  a  heuristic  technique  that  is  robust  to  the 
lack  of  FDS  patches  in  the  ultrasound  image.  This  method  allows  the  calculation  of  the  elevational  distance  for 
all  patches  of  the  image,  regardless  of  their  level  of  coherency,  by  measuring  the  axial  and  lateral  correlation 
of  each  patch.  Since  the  behavior  of  coherent  reflectors  in  the  elevational  direction  can  be  different  from  their 
behavior  in  the  axial  and  lateral  directions,  the  performance  of  the  method  can  decline  depending  on  the  level 
of  anisotropy  of  the  tissue. 

In,14  we  proposed  a  fast  algorithm  to  find  irregularly  shaped  FDS  patches  and  showed  that  this  algorithm 
finds  significantly  more  FDS  patches.  Here,  we  use  beam  steering  as  another  technique  to  increase  the  number 
of  FDS  patches  found  in  the  image.15  This  is  achieved  by  obtaining  more  data  from  a  certain  region  of  tissue, 
hence  reducing  the  size  of  the  analysis  patch.  Having  found  such  small  FDS  patches,  we  further  use  the  steered 
images  for  better  out-of-plane  (elevational)  motion  estimation. 

Coherent  scattering  causes  the  elevational  distance  measurement  from  the  conventional  correlation  algorithms 
to  be  underestimated.11  Thus,  distance  measurement  is  limited  to  the  patches  of  the  ultrasound  image  that  con¬ 
tain  only  FDS.11  To  completely  determine  the  out-of-plane  degrees  of  freedom  between  two  planes,  at  least  three 
non-collinear  pairs  of  such  patches  are  required.4 


Since  FDS  patches  are  extremely  rare  in  real  tissue,  these  methods  usually  have  a  low  accuracy  and  are  only 
relevant  in  limited  tissue  types.  Gee  et  al.13  proposed  a  heuristic  technique  that  is  robust  to  the  lack  of  FDS 
patches  in  the  ultrasound  image.  This  method  allows  the  calculation  of  the  elevational  distance  for  all  patches 
of  the  image,  regardless  of  their  level  of  coherency,  by  measuring  the  axial  and  lateral  correlation  of  each  patch. 
Since  the  behavior  of  coherent  reflectors  in  the  elevational  direction  can  be  different  from  their  behavior  in  the 
axial  and  lateral  directions,  the  performance  of  the  method  can  decline  depending  on  the  level  of  anisotropy  of 
the  tissue.  The  purpose  of  this  work  is  to  devise  a  method  applicable  to  a  various  tissue  types  that  accurately 
reconstructs  3D  volumes  from  ultrasound  images. 

Recently,  Laporte  and  Arbel16, 17  have  proposed  probabilistic  fusion  of  noisy  out-of-plane  motion  estimation. 
This  work  is  most  similar  to  these  works,  in  that  it  calculates  the  maximum  likelihood  estimate  (MLE)  of  the 
out-of-plane  motion.  We  also  use  beam  steering  to  obtain  more  data  and  increase  the  accuracy  of  the  out  of 
plane  motion  estimation  similar  to.15 


2.  METHODS 

2.1.  Combining  Steered  Images 

We  are  looking  for  rectangular  FDS  patches  using  images  acquired  from  the  same  location  at  different  steering 
angles.  The  key  idea  is  to  combine  data  acquired  from  a  certain  region  at  different  steering  angles  and  therefore 
reducing  the  size  of  the  analysis  patch.  Figure  2  shows  two  images  acquired  at  0  and  0  steering  angles.  A 
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Figure  1.  (a)  shows  the  three  directions  relative  to  the  ultrasound  probe.  Out-of-plane  direction  and  elevational  direction 
are  used  interchangably  in  this  work,  (b)  shows  aquisition  of  two  ultrasound  images  at  a  distance  of  Az.  Ultrasound 
beam  is  in  order  of  a  millimeter  wide.  This  wideness  affects  the  resolution  of  ultrasound  image  in  the  lateral,  y,  and 
elevational,  z,  directions,  as  well  as  creating  a  granular  pattern,  called  speckle.  The  size  of  the  resolution  cell  in  the  axial 
direction,  x ,  is  determined  by  the  wavelength  of  the  ultrasound  wave  and  is  magnified  in  this  image. 


rectangle  patch  in  the  left  image  is  warped  into  a  parallelogram  and  is  shifted  in  the  steered  right  image.  The 
position  of  the  parallelogram  can  be  simply  found  as  a  function  of  x  and  y.  Therefore,  samples  nx  and  ny 
from  the  steered  image  correspond  to  samples  nx  and  ny  from  the  non-steered  image  and 


nx 


ny 


Vus  n  .  ,  . 
nx  -  •  —  •  sm(0)  •  ny 

ny 

cos  (0) 


(3) 


where  vus  =  1540000mm/s  is  the  speed  of  ultrasound  in  tissue,  v  is  the  sampling  frequency  of  the  ultrasound 
machine,  n  is  the  total  number  of  the  A-lines  and  w  is  the  width  of  image  in  mm.  To  find  the  correspondence 
of  a  patch,  the  correspondent  of  its  four  corners  are  found  using  these  equations  and  applying  nearest  neighbor 
interpolation.  The  parallelogram  connecting  these  four  corners  is  the  correspondent  of  the  patch. 


2.2.  Maximum  Likelihood  Motion  Estimation 

Assume  we  have  two  parallel  ultrasound  images  with  ground  truth  out-of-plane  distance  z  (Figure  1),  and  that 
we  have  measured  correlation  coefficients  pi  for  i  =  1  •  •  -n  patches  between  the  two  images  (Figure  3).  The  goal 
is  to  find  pz  which  is  maximum  likelihood  estimate  of  2  given  all  the  pi  measurements.  Let  pi  =  fi(zi)  be  the 
calibration  function  that  relates  the  out-of-plane  motion  Zi  to  correlation  coefficient  pi  for  patch  i  (each  patch 
has  a  different  calibration  function  depending  on  its  depth,  see  Figure  6).  Assuming  that  pi  is  drawn  from  a 
Gaussian  distribution  with  mean  f(pz)  and  variance  cr^,  the  conditional  probability  of  pi  is 
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Figure  2.  Corresponding  patches  in  images  acquired  with  different  steering  angles.  In  the  left,  a  patch  is  shown  in  the 
not-steered  image.  In  the  middle,  the  patch  which  corresponds  to  the  same  tissue  is  shown  in  the  scan-converted  steered 
image.  In  the  right,  the  patch  is  shown  in  the  raw  steered  image  (not  scan-converted). 


Figure  3.  The  correlation  coefficient  pi  is  calculated  between  n  patches  of  the  two  images. 


We  assume  that  pi  measurements  are  independent.  Therefore,  the  conditional  probability  of  observing  all  the 
pi  values  will  be  simply  their  multiplication.  Looking  at  the  product  as  a  function  of  p  and  &i  and  taking  its 
logarithm  to  convert  multiplication  to  summation,  we  have  the  familiar  log-likelihood  equation 


L(p\fiz,a2)  =  -^=1 


1  1  2  .  (pi-fi(Pz))2 

2  ^  CTi+ - ^ - 


+  |log(27r) 


(5) 


where  P  and  a2  are  two  vectors  containing  all  the  pi  and  a2  measurements.  In  the  above  equation,  pi  is  the 
correlation  of  two  corresponding  patches  and  is  known,  a2  is  also  known:  it  is  the  variance  of  the  correlation 
and  is  calculated  in  the  calibration  process  (Figure  6).  To  find  the  ML  estimate  of  the  pz,  we  differentiate  this 
equation  with  respect  to  pz  and  set  it  to  zero,  arriving  at 


=0  (6) 

where  f  denotes  derivative  of  /.  Unfortunately  this  equation  is  not  easy  to  solve  for  pz.  Instead,  lets  transform 
pi  to  Zi  and  write  the  log-likelihood  functions  in  terms  of  zi  s.  Equation  5  becomes 


z  |  Pz 
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where  azi  =  jrppp  is  the  transformed  variance  in  the  p  domain  (cr,)  to  the  z  domain  ( azi ).  Differentiating  with 
respect  to  pz  and  setting  it  to  zero  gives 


=  0 


Zi  -  p, 


(8) 


which  can  be  easily  solved  to  give 


Vz  = 


(9) 


Finally,  we  utilize  constraints  in  the  images  to  enhance  out-of-plane  motion  estimation.  Many  surgical 
procedures  such  as  biopsy,  drug  delivery  and  brachytherapy  involve  inserting  a  needle  into  the  tissue.  The 
prior  of  needle  continuity  can  be  used  to  decrease  the  variance  of  the  measured  out  of  plane  motion  (we  are 
assuming  that  the  needle  crosses  US  image  plane  and  is  not  parallel  to  the  image).  Assume  that  the  tip 
of  the  needle  can  be  measured  at  each  image  with  a  variance  of  [cr*eedZJ2,  and  that  the  angle  of  the  needle 
with  the  normal  of  the  ultrasound  image  (i.e  the  angle  between  the  needle  and  the  axes  y  in  Figure  1)  is 
a.  Also,  let  cr^or  denote  the  variance  of  the  out-of-plane  motion  estimation  using  the  correlation  method  and 

a  needle  =  &  needle!  fan(<^)  •  Assuming  both  noises  are  Gaussian,  variance  of  the  final  estimate  which  combines 

2  2 

the  two  estimates  is  a2needle°c°r-  .  It  can  be  easily  shown  that  this  quantity  is  less  than  both  o-^eedle  and  a2cor  , 

a needle  'a cor.  «-  t 

meaning  that  the  resulting  variance  is  less  than  both  initial  variances. 


2.3.  Calibration  and  Data  Acquisition 

The  system  operates  in  two  distinct  modes  -  calibration  mode  and  image-based  3DUS  reconstruction  mode 
(Figure  4).  Both  will  be  described  from  a  process  flow  perspective.  In  the  calibration  mode,  information  necessary 
to  calibrate  the  distance  estimations  is  collected  (Figure  4).  To  this  end,  the  robot  control  component  steps  the 
robot  through  a  series  of  precisely  defined  positions  and  triggers  the  acquisition  of  a  single  US  frame  (RF  data) 
at  each  position  from  a  homogeneous  fully  developed  speckle  (FDS)  phantom.  These  frames  are  associated  with 
their  respective  coordinates  and  stored  for  offline  use.  Then,  the  software  system  reads  the  batch  of  frames  and 
positions  and  subdivides  the  frames  into  distinct  subpatches.  Pairs  of  patches  from  the  same  location  originating 
from  frames  at  different  distances  are  correlated,  thus  creating  a  set  of  (strictly  monotonous)  calibration  (or 
decorrelation)  curves  x,y(d).  These  curves  depend  on  the  characteristics  of  the  selected  probe,  the  imaging 
frequency,  and  the  image  location  x,  y  (in  particular  the  depth  y)  of  the  respective  patches.  Currently,  the  offline 
calibration  process  takes  2-3  minutes  including  scan  time  to  generate  the  needed  calibration  curves  (decorrelation 
curves).  Before  this  recent  development,  manual  data  collection  and  offline  processing  using  MATLAB  scripts 
used  to  take  ^  ofw-r 
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Figure  4.  The  data  acquisition  and  calibration  system  (some  of  the  blocks  are  the  subject  of  future  work). 


(a)  Calibration  experiment 


(b)  Out-of-plane  motion  estimation 
using  real  tissue 


Figure  5.  The  experimental  setup  for  moving  the  probe  out-of-plane  and  acquiring  ultrasound  images.  In  (a),  the  robot 
moves  the  probe  in  the  out-of-plane  direction  while  the  ultrasound  is  imaging  a  FDS  phantom  to  generate  the  calibration 
curves.  In  (b),  the  robot  is  moving  the  probe  while  the  ultrasound  is  imaging  real  tissue,  so  that  the  speckle  correlation 
results  can  be  compared  to  ground  truth  (i.e.  robot  readings). 


3.  EXPERIMENTAL  SETUP 

The  ultrasound  RF  data  was  sampled  with  a  robot-based  system  in  order  to  achieve  reliable,  high-accuracy 
ground  truth  readings  for  the  displacements.  This  will  give  us  the  images  we  need  for  calibration  and  also  the 
gives  us  the  ground  truth  when  we  reconstruct  the  volume.  This  process  yielded  a  series  of  planar-parallel  RF 
slices  through  the  respective  phantom,  in  a  fashion  somewhat  comparable  to  a  freehand  sweep.  The  phantoms 
were  positioned  within  the  workspace  of  a  high-precision  three  degrees-of-freedom  (DoF)  cartesian  robot  stage 
(DMC-21x3  with  three  servo  motor  stages,  by  Galil  Motion  Control;  relative  accuracy  better  0.005  mm).  For 
calibration,  the  stage  translated  the  probe  to  new  positions  every  Ax  =  0.05  mm  apart,  then  triggered  RF 
slice  acquisition  via  a  TTL  signal  connected  to  the  ultrasound  machine’s  ECG  trigger  port,  where  the  data  was 
written  to  file.  For  calibration  data  acquisition,  a  FDS  phantom  is  imaged.  For  volume  reconstruction,  real 
tissue  (beef  steak)  is  used.  Figure  5  shows  the  experimental  setup. 

An  Ultrasonix  ultrasound  machine  (Burnaby,  BC)  with  a  sampling  frequency  ofu  =  20MHz  is  used  to  acquire 
RF  data.  To  calibrate  the  rate  of  image  decorrelation  with  out-of-plane  motion,  RF  data  of  5x80  parallel  frames 
were  acquired  from  a  FDS  phantom  at  an  elevational  distance  of  0.05  mm  between  consecutive  images:  five 
frames  at  each  location  with  steering  angles  of  -5,  -2.5,  0,  2.5  and  5.5  degrees.  The  experimental  setup  is  shown 
in  Figure  5:  the  probe  is  moved  with  a  micrometer  with  the  accuracy  of  .005  mm.  Calibration  results  showed 
that  the  decorrelation  rate  is  not  affected  by  beam  steering. 

Out-of-plane  motion  estimation  was  performed  on  ex- vivo  beef  steak  tissue.  4x80  RF  frames  at  an  elevational 
distance  of  0.05  mm  between  consecutive  frames  were  acquired  using  the  setup  shown  in  Figure  5:  four  images 
at  each  location  with  —5°,  —2.5°,  0°,  2.5°  and  5°  steering  angles. 


Figure  6.  Two  calibration  curves  at  depths  of  10  mm  and  20  mm  and  their  variances.  Note  that  the  calibration  curve 
at  the  deeper  location  drop  slower  with  the  ont-of-nlane  motion. 
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Figure  7.  Subdivision  of  an  image  of  the  FDS  phantom  into  patches  is  shown. 


4.  OUT-OF-PLANE  MOTION  ESTIMATION 

The  correlations  are  calculated  using  Pearsons  linear  correlation  coefficient  p 


p(W1Z)  = 


_ ZwjZi  -  NjMvPz _ 
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where  Wi  and  i  m  1  •  •  •  A,  are  the  intensity  values  of  each  pixel  in  patches  W  and  Z,  N  is  the  total  number 
of  pixels  and  pw  and  py  are  the  means  of  the  intensity  values  of  patches  W  and  Z  respectively. 

Patches  that  are  closest  to  being  FDS  are  selected  as  described  in.15  Figure  8  shows  the  results  of  recon¬ 
structing  out-of-plane  motion  using  the  correlation  values,  (a)  and  (b)  are  obtained  by  combining  the  two  images 
with  ±2.5°  steering  angle  at  each  location,  while  (c)  and  (d)  are  obtained  by  combining  the  two  images  with 
±5°  steering  angle  at  each  location.  The  results  show  that  using  the  MLE  method  slightly  reduces  both  the 
underestimation  error  and  the  variance  of  the  out-of-plane  measurements.  IN  addition,  it  can  be  seen  from  (b) 
and  (d)  that  the  needle  constraint  reduces  the  variance  in  the  measurements. 


5.  DISCUSSION  AND  CONCLUSION 

Out-of-plane  motion  estimation  is  only  studied  here  for  a  fixed  distance  between  two  frames,  0.4mm.  A  study 
of  accuracy  as  the  distance  varies  gives  insight  for  optimum  frame  selection.18, 19  In  freehand  experiments  the 
images  are  not  parallel  as  they  are  in  our  experiments,  and  therefore  the  rotations  between  the  images  need 
to  be  found.11, 13,20  We  showed  before15  that  performing  beam  steering  significantly  increases  the  accuracy  of 
out-of-plane  motion  estimation.  In  this  work,  we  showed  that  MLE  can  also  be  used  to  enhance  the  out-of-plane 
motion  estimation. 
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(a)  relative  error,  ±2.5°  steering  angle 


(b)  standard  deviation,  ±2.5°  steering  angle 
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Figure  8.  relative  error  and  standard  deviation  of  the  sensorless  measurements,  (a)  The  relative  error.  Reconstruction 
is  performed  using  two  steered  images  at  ±2.5°.  (b)  The  standard  deviation  of  the  measurements.  W/O  MLE  refers  to 
without  MLE,  W/O  constraint  refers  to  without  utilizing  the  needle  continuity  constraint,  and  W  MLE  refers  to  with 
MLE.  (c)  and  (d)  are  corresponding  errors  and  variances  with  two  steered  images  at  ±5°. 
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Abstract.  Tracked  ultrasound  elastography  can  be  used  for  guidance  in 
partial  breast  radiotherapy  by  visualizing  the  hard  scar  tissue  around  the 
lumpectomy  cavity.  For  clinical  success,  the  elastography  method  needs 
to  be  robust  to  the  sources  of  decorrelation  between  ultrasound  images, 
specifically  fluid  motions  inside  the  cavity,  change  of  the  appearance  of 
speckles  caused  by  compression  or  physiologic  motions,  and  out-of-plane 
motion  of  the  probe.  In  this  paper,  we  present  a  novel  elastography  tech¬ 
nique  that  is  based  on  analytic  minimization  of  a  regularized  cost  func¬ 
tion.  The  cost  function  incorporates  similarity  of  RF  data  intensity  and 
displacement  continuity,  making  the  method  robust  to  small  decorre¬ 
lations  present  throughout  the  image.  We  also  exploit  techniques  from 
robust  statistics  to  make  the  method  resistant  to  large  decorrelations 
caused  by  sources  such  as  fluid  motion.  The  analytic  displacement  esti¬ 
mation  works  in  real-time.  Moreover,  the  tracked  data,  used  for  targeting 
the  radiotherapy,  is  exploited  for  discarding  frames  with  excessive  out- 
of-plane  motion.  Simulation,  phantom  and  patient  results  are  presented. 


1  Introduction 

Breast  irradiation  after  lumpectomy  significantly  reduces  the  risk  of  cancer  re¬ 
currence.  There  is  growing  evidence  suggesting  that  irradiation  of  only  the  in¬ 
volved  area  of  the  breast,  partial  breast  irradiation  (PBI),  is  as  effective  as  whole 
breast  irradiation  [1].  Benefits  of  PBI  include  significantly  shortened  treatment 
time  and  fewer  side  effects  as  less  tissue  is  treated.  However,  these  benefits  cannot 
be  realized  without  localization  of  the  lumpectomy  cavity.  Tracked  ultrasound 
elastography  can  be  used  for  localizing  the  lumpectomy  cavity  in  the  treatment 
room,  minimizing  tissue  motion  from  planning  to  treatment. 

This  paper  is  focused  on  freehand  palpation  elastography,  which  involves  es¬ 
timating  the  displacement  field  of  the  tissue  undergoing  slow  compression.  Most 
elastography  techniques  estimate  the  displacement  field  using  local  cross  corre¬ 
lation  analysis  of  echoes  [2,3,4].  These  methods  are  very  sensitive  and  accurate 
for  calculating  small  displacements.  However,  elastography  is  subject  to  speckle 
decorrelation  caused  by  various  sources  such  as  motion  of  subresolution  scatter¬ 
ed,  out-of-plane  motion,  high  compression  and  complex  fluid  motions. 

The  prior  of  tissue  deformation  continuity  can  be  used  to  make  elastography 
more  robust  to  signal  decorrelation.  Previous  work  on  regularized  elastography 
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is  computationally  expensive  [5,6].  Dynamic  programming  (DP)  can  be  used 
to  speed  the  optimization  procedure  [7],  but  it  only  gives  integer  displacements. 
Subpixel  displacement  estimation  is  possible  [7],  but  it  is  computationally  expen¬ 
sive  if  a  fine  subpixel  level  is  desired.  In  addition,  a  fixed  regularization  weight  is 
applied  throughout  the  image.  However,  while  two  ultrasound  images  may  cor¬ 
relate  well  in  most  parts,  they  can  have  small  correlation  in  specific  parts.  Four 
examples  of  low  correlation  are:  (1)  correlation  decreases  with  depth  mainly  due 
to  a  decrease  in  the  ultrasonic  signal  to  noise  ratio,  (2)  correlation  is  low  close 
to  arteries  due  to  complex  motion  and  inside  vessels  due  to  blood  motion,  (3) 
correlation  is  extremely  low  in  lesions  that  contain  liquid  due  to  the  incoherent 
fluid  motion  [8,3],  and  (4)  out-of-plane  motion  of  movable  structures  within  the 
image  [8]  causes  low  local  correlation.  To  prevent  such  regions  from  introduc¬ 
ing  errors  in  the  displacement  estimation  one  should  use  large  weights  for  the 
regularization  term,  resulting  in  over-smoothing. 

Freehand  palpation  elastography  provides  ease-of-use  and  requires  minimum 
additional  cost.  However,  out-of-plane  motion  cannot  be  avoided  in  freehand  pal¬ 
pation,  which  reduces  the  quality  of  any  elastography  method.  Assisted  freehand 
elastography  [9]  significantly  reduces  the  out-of-plane  motion  but  it  requires  ad¬ 
dition  of  a  device  to  the  probe.  Quality  metrics  such  as  persistence  in  strain 
images  have  also  been  developed  to  address  this  problem  [10].  To  measure  the 
persistence,  elastography  is  performed  on  two  pairs  of  images  and  the  resulting 
strain  images  are  correlated.  This  method  requires  strain  images  for  calculating 
the  quality  metric.  Therefore,  trying  all  the  combinations  in  a  series  of  frames 
to  find  the  best  pair  for  elastography  will  be  computationally  expensive. 

In  this  paper,  we  present  a  novel  elastography  method  based  on  analytic 
minimization  (AM)  of  a  cost  function  that  incorporates  similarity  of  echo  ampli¬ 
tudes  and  displacement  continuity.  We  introduce  a  novel  regularization  term  and 
demonstrate  that  it  minimizes  displacement  underestimation  caused  by  smooth¬ 
ness  constraint.  We  also  introduce  the  use  of  robust  statistics  implemented  via 
iterated  reweighted  least  squares  (IRLS)  to  treat  uncorrelated  ultrasound  data 
as  outliers.  And  finally,  we  use  the  tracking  information  to  select  the  best  pairs 
of  frames  for  elastography.  Simulation,  phantom  and  patient  experiments  are 
presented  for  validation. 

2  Regularized  Displacement  Estimation 

Dynamic  Programming  (DP).  DP  is  a  discrete  efficient  optimization  tech¬ 
nique  for  causal  systems.  In  DP  elastography  [7],  a  cost  function  is  defined  as 

C(i,di)  =  min  {C(i  -  l,d*-i)  +  aaR(di,  di-i)}  +  \h(i)  -  I2(i  +  di) |  ,  i  =  2---m 

di- 1 

(1) 

where  di  is  the  displacement  of  sample  i,  R{di,di-\)  =  (di  —  di- 1)  is  an  axial 
regularization  term  (axial,  lateral  and  out-of-plane  directions  are  respectively 
z,  x  and  y  in  Figure  2  (a)),  aa  is  a  weight  for  the  regularization,  I\  and  I2 
are  corresponding  RF-lines  of  before  and  after  deformation  and  m  is  the  length 
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of  RF-lines.  The  cost  function  is  minimized  at  i  =  m  and  the  di  values  that 
have  minimized  the  cost  function  are  traced  back  to  i  =  1,  giving  the  di  for  all 
samples.  We  have  implemented  a  2D  DP  algorithm  similar  to  [7]  to  generate 
integer  displacements  as  a  starting  point  for  the  next  step  of  our  algorithm. 

Analytic  Minimization  (AM).  We  now  propose  a  method  that  analytically 
minimizes  a  regularized  cost  function  and  gives  the  refined  displacement  field. 
Only  axial  displacements  will  be  refined  for  strain  calculation. 

Having  the  integer  displacements  di  from  DP,  it  is  desired  to  find  Adi  values 
such  that  di  +  Adi  gives  the  value  of  the  displacement  at  the  sample  i  for  i  = 
1  •  •  •  m.  Such  Adi  values  will  minimize  the  following  regularized  cost  function 

C  (Adi,  •  •  • ,  Adm )  =  r™!  [h(i)  -I2(i  +  di  +  Adi)}2  + 
aa(di  +  Adi  ~  di- $  -  Adi- 1)2  +  a/(d»  +  Adt  -  df  -  Ad\')2  (2) 

where  superscript  p.  refers  to  the  previous  RF-line  (adjacent  RF-line  in  the 
lateral  direction)  and  ai  is  a  weight  for  lateral  regularization.  Substituting  l2(i  + 
di~\~  Adi)  with  its  first  order  Taylor  expansion  approximation  around  di,  we  have 

C  {Adi,  •  •  • ,  Adm)  =  £™1  lh(i)  -  I2(i  +  di)  -  /'(*  +  di)Adi )f  + 
aa(di  +  Adi  ~  di-i  -  Adi- i)2  +  «/(<**  +  Adt  ~  df  -  Adf)2  (3) 

where  J2  is  the  derivative  of  the  I2 .  The  optimal  Adi  values  occur  when  the 
partial  derivative  of  C  w.r.t.  Adi  is  zero.  Setting  dd^d  =  0  we  have 


(I'22  +  aaT>  +  ati)Ad  =  I'2e  -  (aaD  +  azi)d  +  ,  D  = 


1  -1  0  •••  0 

-1  2  -1  •••  0 


[0  •••  0  -1  1J 

(4) 

where  I2  =  diag(/2(l  +  di)  •  •  •  J2(ra  +  dm)),  Ad=[Adi  •  •  •  Adm]  ,  e=  [e\  •  •  •  em]  , 
ei  =  I\(i)  —  12(1  +  di),  d  =  [di  *  •  •  dm]T,  d t-p-  =  dp  +  Adp  is  the  vector  of  total 
displacement  of  the  previous  line  and  I  is  the  identity  matrix.  I2,  D  and  I  are 
matrices  of  size  m  x  m  and  Ad ,  r,  d  and  dt  p-  are  vectors  of  size  m. 


Biasing  the  Regularization.  The  regularization  term  aa(di  +  Adi  —  di-i  ~ 
Adi- 1)2  penalizes  the  difference  between  di~\~Adi  and  di-\+Adi-\,  and  therefore 
can  result  in  underestimation  of  the  displacement  field.  Such  underestimation 
can  be  prevented  by  biasing  the  regularization  by  e  to  aa(di  +  Adi  ~  di-i  — 
Adi- 1  —  e)2,  where  e  =  (dm  —  di)/(rrt  —  1)  is  the  average  displacement  difference 
between  samples  i  and  i  —  1.  An  accurate  enough  estimate  of  dm  —  d\  is  known 
from  the  previous  line.  With  the  bias  term,  the  R.H.S.  of  Equation  4  becomes 
I2e  —  (<aaD  +  ail)d  +  o/dt  ;p-  +  b  where  the  bias  term  is  b  =  aa[— e  0  •  •  •  0  e]T 
and  all  other  terms  are  as  before.  Interestingly,  except  for  the  first  and  the  last 
equation  in  this  system,  all  other  m  —  2  equations  are  same  as  Equation  4. 

Equation  4  can  be  solved  for  Ad  in  4m  operations  since  the  coefficient  matrix 
I22  +  aa~D  -\-aiI  is  tridiagonal.  Utilizing  its  symmetry,  the  number  of  operations 


510 


H.  Rivaz  et  al. 


can  be  reduced  to  2m.  The  number  of  operations  required  for  solving  a  system 
with  a  full  coefficient  matrix  is  more  than  m3/ 3,  significantly  more  than  2m. 

Making  Tracking  Resistant  to  Outliers.  Even  with  pure  axial  compression, 
some  regions  of  the  image  may  move  out  of  the  imaging  plane  and  increase 
the  decorrelation.  In  such  parts  the  confidence  of  the  data  term  is  less  and 
therefore  the  weight  of  the  regularization  term  should  be  increased.  The  parts 
of  the  image  with  low  correlation  can  be  regarded  as  outliers  and  therefore 
a  robust  estimation  technique  can  limit  their  effect.  Before  deriving  a  robust 
estimator  for  Z\d,  we  rewrite  Equation  3  as  C(Ad)  =  T’]T1p(r^)  +  R(Ad)  where 
ri  =Ii{i)  —  l2{iJrdi)  —  H2{iJtdi)Adi1p{ri)  =  rf  and  R  is  the  regularization  term. 
The  M-estimate  of  Ad  is  Ad  =  arg ruined  (T’]T1p(r^)  +  R(Ad)}  where  p(u)  is 
a  robust  loss  function  [11].  The  minimization  is  solved  by  setting  dd^d.  =  0: 


p\n) 


dr 

dAdi 


dR(Ad) 

dAdi 


(5) 


A  common  next  step  [11]  is  to  introduce  a  weight  function  re,  where  w(ri).ri  = 
p'(ri).  This  leads  to  a  process  known  as  “iteratively  reweighted  least  squares” 
(IRLS),  which  alternates  steps  of  calculating  weights  w(ri)  for  r*  =  1  •  •  •  m  using 
the  current  estimate  of  Ad  and  solving  Equation  5  to  estimate  a  new  Ad  with 
the  weights  fixed.  Among  many  proposed  shapes  for  w(-),  we  use  [11] 


w(ri) 


i  \n\<T 

w\  N>T 


(6) 


where  T  is  a  threshold  that  can  be  tuned.  A  small  T  will  treat  many  samples  as 
outliers.  With  the  addition  of  the  weight  function,  Equation  5  becomes 


(wl^2  +  aD  +  a2l)Ad  =  wl^e  —  (oqD  +  o^ijd  +  a2dt'p'  +  b  (7) 


where  w  =  diag(w(r\)  •  •  -w(rm)).  All  of  the  results  presented  in  this  work  are 
obtained  with  one  iteration  of  the  above  equation  unless  otherwise  specified. 
Current  implementation  of  the  AM  algorithm  with  the  IRLS  takes  0.015s  to 
generate  a  dense  displacement  field  of  size  1300  x  60  on  a  3.4GHz  P4  CPU(not 
including  the  DP  run  time).  The  computation  time  increases  linearly  with  the 
size  of  images. 

Frame  Selection.  The  ultrasound  probe  is  tracked  in  navigation/guidance 
systems  to  provide  spatial  information,  to  generate  freehand  3D  ultrasound, 
or  to  facilitate  multi-modality  registration.  Through  a  calibration  process,  the 
6DOF  motion  of  the  probe  in  the  sensor  coordinate  system  is  transformed  into 
image  coordinate  system  [12].  The  mean  of  the  absolute  motion  value  of  all  pixels 
in  3D,  (1^1),  {\vy\)  and  (\vz\),  can  be  analytically  related  to  the  6DOF  sensor 
readings  using  straightforward  and  efficient  geometric  computations.  For  frame 
i  and  j  to  be  selected  from  a  sequence  of  frames  for  elastography, 


Qi 


kX  (I'Lrl)2  +  ky  ( \vy\ )2  +  kz 


<N> 


Uz,opt  | 


(N>  +  c 


(8) 
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should  be  minimized  where  kx,  ky,  and  kz  are  weights  for  lateral,  out-of-plane 
and  axial  displacements  and  vz^opt  is  the  optimum  axial  motion.  Please  refer  to 
[12]  for  a  rationale  of  the  shape  this  function.  Note  that  the  selected  pairs  are 
not  necessarily  consecutive  frames.  The  parameters,  kx,  ky,  kz ,  vZj0pt  and  c  are 
manually  tuned  to  1,  2,  1,  0.7  and  1  for  the  AM  elastography  method. 

3  Simulation,  Phantom  and  Patient  Results 

Simulation  Results.  RF  ultrasound  data  of  two  phantoms  are  simulated  using 
Field  II  [13].  The  first  phantom  is  50  x  10  x  55mm  and  the  second  one  is  36  x 
10  x  25mm.  They  are  both  made  of  homogeneous  and  isotropic  material:  the 
first  one  is  uniform  and  the  second  one  contains  a  circular  hole  filled  with  water, 
simulating  a  blood  vessel  in  tissue  (Figure  2  (a)).  A  uniform  compression  in  the 
z  direction  is  applied  and  the  3D  displacement  field  of  the  phantom  is  calculated 
using  ABAQUS  finite  element  package  (Providence,  RI).  The  Poisson’s  ratio  is 
set  to  v  =  0.49  in  both  phantoms  to  mimic  real  tissue,  which  causes  the  phantom 
to  deform  in  x  &  y  directions  as  a  result  of  the  compression  in  the  z  direction. 

Respectively  5  x  105  and  1.4  x  105  scatterers  with  uniform  scattering  strengths 
are  uniformly  distributed  in  the  first  and  second  phantom,  ensuring  more  than 
10  scatterers  exist  in  a  resolution  cell.  The  scatterers  are  distributed  in  the  8mm 
diameter  vein  also  (Figure  2  (a)).  To  construct  deformed  ultrasound  images,  the 
displacement  of  all  of  the  scatterers  is  calculated  by  interpolating  the  displace¬ 
ment  of  the  neighboring  nodes  in  the  finite  element  analysis.  The  parameters 
of  the  probe  are  set  to  mimic  Siemens  5-10MHz  probes.  The  probe  frequency  is 
7.27MHz,  the  sampling  rate  is  40MHz  and  the  fractional  bandwidth  is  60%. 

The  first  phantom  undergoes  uniform  compressions  in  the  z  direction  to 
achieve  strain  levels  of  2%  to  14%  in  2%  intervals.  Ground  truth  integer  dis¬ 
placement  values  are  used  as  the  initial  estimate  for  AM  to  decouple  the  perfor¬ 
mance  of  DP  from  AM.  Accurate  subpixel  displacement  field  is  calculated  with 
AM  and  the  mean  strain  values  are  compared  with  the  ground  truth  (Figure  1 
(a)-(c)).  The  results  are  only  shown  for  2%,  4%,  8%  and  14%  compression  for 
better  visualization.  The  results  with  two  threshold  values  for  IRLS  and  without 
IRLS  demonstrate  that  outlier  rejection  does  not  affect  the  mean  strain  value, 
while  increasing  the  regularization  weight  aa  increases  underestimation  of  the 
displacement.  The  rate  of  increase  of  the  underestimation  with  increasing  aa  is 
significantly  more  with  the  unbiased  regularization  (dashed  line)  as  expected. 

Significantly  higher  signal  to  noise  ratio  (SNR)  [2]  values  can  be  achieved 
with  outlier  rejection  (Figure  1  (d)-(f))  without  over-smoothing  the  image  with 
high  aa  values.  To  show  the  performance  of  the  overall  method,  the  initial  inte¬ 
ger  displacement  field  is  calculated  with  DP  and  accurate  displacement  field  is 
calculated  with  (Figure  1  (g)-(i)).  The  SNR  values  are  less  than  previous  case 
especially  at  high  strain  values,  where  DP  results  deviates  from  ground  truth. 

The  second  simulation  experiment  is  designed  to  show  the  effect  of  smooth¬ 
ness  weight  and  IRLS  threshold  on  contrast  to  noise  ratio  (CNR)  [2]  when  the 
correlation  is  lower  in  parts  of  the  image  due  to  fluid  motion.  The  phantom 
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a  a  a 

a  a  a 


(a)  T  =  0.005  (b)  T  =  0.01  (c)  no  IRLS 


(d)  T  =  0.005  (e)  T  =  0.01  (f)  no  IRLS 


(g)  T  =  0.005  (h)  T  =  0.01  (i)  no  IRLS 

Fig.  1.  Mean  and  SNR  of  the  elastograms  of  the  Field  II  simulated  uniform  phantom 
at  four  different  compression  levels  (shown  in  percentage)  for  three  IRLS  T  values.  The 
solid  and  dashed  lines  correspond  to  biased  and  unbiased  regularizations  respectively. 
(a)-(c)  shows  the  relative  underestimation  of  the  strain,  e  is  the  mean  strain  calculated 
with  the  elastography  method  and  e*  is  the  ground  truth,  (d)-(f)  shows  the  SNR  of 
the  AM.  (g)-(i)  shows  the  SNR  of  the  AM  with  initial  displacements  found  by  DP. 


contains  a  vein  oriented  perpendicular  to  the  image  plane  (Figure  2).  The  initial 
integer  displacement  is  generated  with  DP.  The  background  window  for  CNR 
calculation  is  located  close  to  the  target  window  to  show  how  fast  the  strain  is 
allowed  to  vary,  a  property  related  to  the  spatial  resolution.  The  maximum  CNR 
with  IRLS  is  5.3  generated  at  T  =  0.005  and  aa  =  38,  and  without  IRLS  is  4.8 
at  aa  =  338.  Such  high  aa  value  makes  the  share  of  the  data  term  in  the  cost 
function  very  small  and  causes  over-smoothing. 

Phantom  Results.  We  perform  freehand  palpation  experiment  on  a  breast 
phantom  to  examine  the  performance  of  the  frame  selection  technique.  50  frames 
of  RF  ultrasound  data  are  acquired  using  a  Siemens  Antares  system  (Issaquah, 
WA).  Our  custom  data  acquisition  program  is  connected  to  the  Axius  Direct 
Research  Interface  to  send  the  command  for  capturing  RF  data.  At  the  same 
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(a)  simulation  phan-  (b)  finite  element  strain 
tom 


(c)  CNR 


Fig.  2.  The  target  (circular)  and  background  (rectangular)  windows  for  CNR  calcula¬ 
tion  of  (c)  are  shown  in  (b) 


Frame  number  Frame  number 

Fig.  3.  The  SNR  and  CNR  of  the  phantom  experiment  with  and  without  frame  selec¬ 
tion 


(a)  CT  image 


(b)  B-mode  image 


width  (mm) 

(c)  strain  image 


Fig.  4.  Patient  experiment  results.  The  arrow  points  to  the  lumpectomy  cavity. 


time,  the  program  collects  tracking  information  from  a  Polaris  tracker  (Waterloo, 
Canada) .  Currently,  the  RF  frames  are  stored  on  the  ultrasound  system  and  are 
processed  offline.  Figure  3  shows  the  SNR  and  CNR  results.  In  automatic  frame 
selection,  Qij  (equation  8)  for  any  two  frames  in  a  buffer  of  size  15  frames 
is  calculated.  For  the  two  frames  which  give  the  minimum  Q,  the  strain  image 
is  obtained.  The  next  image  is  then  fed  to  the  buffer,  its  first  image  is  removed 
and  the  frame  selection  is  performed  again.  The  automatic  frame  selection  gives 
8  frame  pairs  for  strain  calculation  (as  seen  in  the  figure  by  8  SNR  and  CNR 
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values).  Without  frame  selection,  49  strain  images  are  calculated.  The  average 
CNR  and  SNR  values  are  improved  from  4.91  to  7.19  and  from  4.98  to  5.88  with 
frame  selection. 

Patient  Results.  We  have  acquired  freehand  palpation  ultrasound  RF  data 
using  the  Siemens  Antares  system  from  patients  approximately  four  weeks  after 
lumpectomy.  The  ultrasound  probe  is  tracked  with  the  Polaris  tracking  system. 
Optimal  frame  selection  is  performed  to  select  images  for  elastography  using  the 
AM  method.  The  strain  image  (Figure  4)  shows  that  the  AM  method  can  detect 
the  thin  hard  scar  tissue  even  though  it  is  close  to  the  cavity  fluids  which  undergo 
incoherent  motions  and  cause  signal  decorrelation.  Since  the  AM  method  finds 
the  displacement  of  all  the  samples  on  an  A-line  at  the  same  time,  the  correlated 
data  at  the  top  and  bottom  of  the  cavity  guide  the  method  to  find  the  correct 
displacement  inside  the  cavity  where  the  data  is  decorrelated. 

4  Discussion  and  Conclusion 

We  introduced  a  novel  method  for  calculating  a  dense  displacement  map  by 
analytic  minimization  of  a  cost  function.  We  used  the  IRLS  method  from  ro¬ 
bust  statistics  to  make  the  tracking  resistant  to  outliers.  Moreover,  we  exploited 
the  tracking  data  to  optimize  frame  selection.  Through  simulation  studies  using 
Field  II  and  finite  element  analysis,  we  showed  that  the  proposed  AM  method 
generates  high  quality  displacement  estimates.  The  elastography  method  works 
in  real-time.  A  comparison  of  the  IRLS  method  with  quality  guided  displacement 
tracking  [14]  which  also  aims  for  robustness  is  a  subject  of  future  work. 

We  chose  the  novel  application  of  the  lumpectomy  cavity  localization  as  the 
hard  scar  tissue  is  relatively  thin  and  demands  a  high  resolution  elastography 
method.  Also,  incoherent  fluid  motions  in  the  cavity  causes  large  decorrelations, 
requiring  a  robust  method.  We  have  an  active  Institutional  Review  Board(IRB) 
protocol  and  have  promising  results  from  9  patients  which  will  be  published  in 
future  work. 
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Abstract.  This  paper  presents  a  robust  framework  for  freehand  ultra¬ 
sound  elastography  to  cope  with  uncertainties  of  freehand  palpation  us¬ 
ing  the  information  from  an  external  tracker.  In  order  to  improve  the 
quality  of  the  elasticity  images,  the  proposed  method  selects  a  few  im¬ 
age  pairs  such  that  in  each  pair  the  lateral  and  out-of-plane  motions  are 
minimized.  It  controls  the  strain  rate  by  choosing  the  axial  motion  to  be 
close  to  a  given  optimum  value.  The  tracking  data  also  enables  fusing 
multiple  strain  images  that  are  taken  roughly  from  the  same  location. 
This  method  can  be  adopted  for  various  trackers  and  strain  estimation 
algorithms.  In  this  work,  we  show  the  results  for  two  tracking  systems 
of  electromagnetic  (EM)  and  optical  tracker.  Using  phantom  and  ex-vivo 
animal  experiments,  we  show  that  the  proposed  techniques  significantly 
improve  the  elasticity  images  and  reduce  the  dependency  to  the  hand 
motion  of  user. 

Key  words:  Ultrasound,  Elastography,  Elasticity,  Tracking,  Strain 


1  Introduction 

Ultrasound  elastography  is  an  emerging  medical  imaging  modality  which  in¬ 
volves  imaging  the  mechanical  properties  of  tissue  and  has  numerous  clinical 
applications.  Among  many  variations  of  ultrasound  elastography  [1],  our  work 
focuses  on  real-time  static  elastography,  a  well-known  technique  that  applies 
quasi-static  compression  of  tissue  and  simultaneously  images  it  with  ultrasound. 
Within  many  techniques  proposed  for  static  elastography,  we  focus  on  freehand 
palpation  elasticity  imaging  which  involves  deforming  the  tissue  by  simply  press¬ 
ing  the  ultrasound  probe  against  it.  Freehand  ultrasound  elastography  has  shown 
great  potential  in  clinical  applications  especially  for  diagnosis  and  screening  of 
breast  lesions  [2].  The  application  of  elastography  is  not  limited  to  breast,  and 
other  applications  such  as  diagnosis  of  prostate  cancer,  monitoring  ablation  and 
deep  vein  thrombosis  have  also  been  studied. 

Despite  the  reports  on  success  of  elastography,  yet  it  has  not  become  a  part  of 
any  routine  clinical  application.  The  main  reason  is  that  elastography  is  highly 
qualitative  and  user-dependent.  The  best  result  is  achieved  when  the  user  com¬ 
presses  and  decompresses  the  tissue  uniformly  in  the  axial  direction  with  the 
proper  hand  motion.  It  is  difficult  to  control  the  compression  rate  as  it  is  gov¬ 
erned  by  the  hand  motion  and  the  frame  rate  of  RF  data.  Also,  small  lateral 
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or  out-of-plane  motions  can  compromise  the  quality  of  images.  However,  it  is 
difficult  to  induce  pure  axial  motion  with  freehand  compression.  Sophisticated 
algorithms  can  only  partially  address  the  problem  by  compensating  for  in-plane 
motions  and  applying  smoothness  constraints.  The  images  are  also  hard  to  in¬ 
terpret,  and  artifacts  -caused  by  failure  of  the  strain  estimation  algorithm  or 
poor  hand  motion-  may  be  mistaken  for  lesions  inside  the  soft  tissue.  Devel¬ 
oping  an  elastography  technique  that  is  not  affected  by  poor  hand  motion  and 
other  sources  of  signal  decorrelation  will  pave  the  way  for  wide-spread  clinical 
use  of  elastography. 

To  improve  the  reliability,  quality  metrics  such  as  persistence  in  strain  images 
have  been  developed  [3,4].  This  quality  indicator  is  calculated  for  each  image 
and  provided  to  the  user  as  feedback.  Persistence  is  also  used  to  merge  multiple 
elasticity  images  together  [3].  To  measure  the  persistence,  strain  is  computed 
for  two  pairs  of  echo  frames,  and  the  resulting  images  are  correlated.  Although 
these  techniques  offers  a  major  advantage,  there  remains  several  limitations. 
First,  the  strain  has  to  be  estimated  before  the  calculation  of  the  quality  metric. 
With  typical  ultrasound  settings,  the  frame  rate  can  reach  more  than  30  Hz.  For 
subsequent  frames,  an  efficient  implementation  of  this  image-based  metric  might 
cope  with  this  rate.  Nonetheless,  the  task  will  be  extremely  difficult  to  try  all 
the  combinations  in  a  series  of  frames.  Moreover,  the  quality  metric  will  not  be 
able  to  provide  feedback  to  the  user  whether  he/she  should  adjust  the  palpation 
in  certain  direction.  Also,  there  would  be  minimal  control  over  the  strain  rate. 

The  ultrasound  probe  is  often  tracked  in  navigation/guidance  systems  to 
provide  spatial  information,  to  form  freehand  3D  ultrasound,  or  to  facilitate 
multi- modality  registration.  In  this  work,  we  exploit  the  tracking  data  to  en¬ 
hance  the  quality  of  the  elasticity  images.  We  use  the  tracking  data  to  select 
multiple  image  pairs  that  contain  the  optimum  deformation  for  the  elastogra¬ 
phy  algorithm.  The  optimum  value  for  lateral  and  out-of-plane  motions  is  zero, 
and  the  optimum  axial  motion  is  determined  by  the  specific  elastography  algo¬ 
rithm  used,  which  is  Normalized  Cross- Correlation  (NCC)  in  this  work.  Next, 
we  fuse  the  strain  images  obtained  from  the  multiple  image  pairs  together  based 
on  the  location  of  each  strain  image  to  improve  image  quality.  We  assume  that 
the  ultrasound  data  is  2D.  Nonetheless  similar  techniques  proposed  here  could 
be  extended  to  3D  ultrasound. 

2  Methodology 

Consider  a  sequence  of  RF  data  collected  during  the  palpation  of  tissue  using  a 
tracked  transducer.  We  have  previously  shown  that  it  is  possible  to  synchronize 
the  RF  frames  with  the  tracking  information  relying  only  on  the  same  data 
collected  during  palpation  [5].  From  synchronization,  the  tracking  information 
is  interpolated  at  the  incident  time  of  each  frame.  The  input  to  our  algorithm  is 
then  a  series  of  RF  frames  along  with  their  corresponding  transformation. 

First,  we  need  to  define  a  distance  function  between  two  frames  of  RF  data. 
For  this  purpose,  we  use  a  model  of  image  decorrelation  in  presence  of  out-of- 
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plane  and  lateral  motion.  RF  signal  is  often  modeled  as  the  collective  response 
of  scatterers  randomly  distributed  within  the  resolution  cell  of  the  ultrasound 
[6,7].  Each  scatterer  is  assumed  to  have  an  amplitude  governed  by  the  shape 
of  the  resolution  cell  and  a  phase  which  is  distributed  from  0  to  tt  uniformly  at 
random.  Considering  a  Gaussian  shape  for  the  resolution  cell  Prager  et.  al  [8] 
calculated  the  correlation  as  a  function  of  out-of-plane  motion  to  be  exp{—^ 2). 
S  and  a  denote  the  displacement  and  the  width  of  the  resolution  cell  respectively. 
Although  this  function  is  only  valid  for  fully  developed  speckle,  it  provides  a  con¬ 
venient  estimate  of  correlation.  It  should  be  noted  that  in  [8] ,  the  displacement  is 
estimated  from  correlation,  whereas  here,  we  intend  to  define  an  energy  function 
based  on  displacement.  Extending  this  formula  to  both  out-of-plane  and  lateral 
displacements,  we  define  our  energy  function,  E(x,z),  as  follows: 


E(DX,  Dz)  =  exp{—Kx  ■  D2X  -  Kz  ■  D2Z),  (1) 

where  Dx  and  Dz  represent  the  displacement  in  out-of-plane  and  lateral  direc¬ 
tions  ( Dy  is  reserved  for  axial  motion).  Kx  and  Kz  determine  the  sensitivity  to  a 
certain  direction.  In  order  to  be  able  to  use  this  function,  we  need  a  component¬ 
wise  metric  representing  the  distance  of  two  frames  given  their  homogeneous 
transformations.  The  first  step  is  to  compute  the  relative  transformation  be¬ 
tween  them.  Suppose  a  =  [ax  ay  az]T  is  the  axis-angle  representation  of  the 
relative  rotation,  and  t  =  [tx  tytz]T  is  the  relative  translation.  Assuming  a  small 
rotation,  the  relative  displacement  of  a  point,  P  =  [xy  0]T,  will  be  d  =  a  x  P  +  t. 
We  then  define  the  distance  vector  of  two  frames,  D  =  [Dx  Dy  DZ]T ,  as  the  RMS 
of  the  components  of  d  for  all  the  points  in  the  region  of  interest  (ROI): 
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where  sqrt{.}  returns  the  root.  Here,  ROI  is  assumed  to  be  rectangular  and 
determined  by  xi,  #2,  2/i,  and  y^.  The  vector  D  provides  a  measure  of  distance 
for  each  direction  separately.  We  use  this  vector  in  Equation  (1)  which  gives  us 
an  estimate  of  “pseudo-correlation”  over  the  ROI. 

The  data  goes  through  four  stages  of  processing  to  create  a  single  high- 
quality  strain  image.  In  the  first  step,  few  images  are  selected  from  the  data 
series  that  are  approximately  collected  from  one  cross-section  of  tissue  with 
minimal  lateral  and  out-of-plane  motion.  To  this  end,  the  energy  function  of 
each  frame  is  computed  with  respect  to  all  other  frames  in  the  sequence.  Then, 
the  total  energy  is  found  for  each  frame  as  the  sum  of  the  energies  of  the  M 
closest  frames,  where  closeness  implies  higher  energy,  and  M  is  the  maximum 
number  of  frames  to  be  selected.  Then,  the  frame  with  the  highest  total  energy 
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(the  center  frame)  is  identified,  and  the  M  closest  frames  to  the  center  frame 
including  itself  are  selected.  Additionally,  the  frames  that  have  E  of  less  than 
0.5  with  respect  to  the  center  frame  are  disqualified.  This  is  applied  to  ensure 
lower  number  of  frames  are  chosen  when  M  frames  from  one  cross-section  are 
not  available. 

In  the  next  stage,  the  program  evaluates  all  possible  combination  of  frame 
pairs  for  elastography.  For  M  frames,  there  will  be  (^)  =  M(M  —  l)/2  pair 
combinations  which  will  be  compared  using  a  slightly  modified  version  of  E. 
Since  the  pairs  are  directly  compared,  it  suffices  to  minimize  the  exponent  of 
Equation  (1)  in  order  to  maximize  E.  We  also  add  a  term  for  axial  motion  that 
penalizes  compressions  that  are  higher  than  an  optimum  compression  value,  topt. 
Hence,  a  “cost  function”,  C 1,  is  defined  as  follows: 

Cl(D)  =  K,  Dl+KrDy2+K,-Dl  j£»  ;  °  (3) 

topt  implies  the  optimal  strain.  Optimal  strain  can  be  theoretically  defined 
as  described  in  [9].  It  also  depends  on  the  robustness  of  the  elasticity  estima¬ 
tion  algorithm.  topt  might  be  within  the  range  of  the  resolution  of  the  tracker. 
Therefore,  at  this  stage  we  do  not  assign  a  penalty  for  the  compressions  less  than 
topt.  If  the  compression  is  close  to  zero,  the  contrast  of  the  reconstructed  image 
degrades.  The  program  filters  the  pairs  with  low  compression  in  the  next  stage 
using  image  content.  Similar  to  the  first  part,  a  maximum  number  of  frames 
with  lowest  cost  are  selected  provided  that  the  cost  is  lower  than  a  threshold. 
The  threshold  is  not  strict  to  ensure  acceptable  pairs  are  not  filtered. 

The  final  pairs  are  selected  by  recovering  the  global  lateral  motion  and  com¬ 
pression  by  matching  the  two  RF  frames  in  each  pair.  The  tracking  information 
is  used  to  initialize  the  search.  For  instance,  the  search  range  for  compression  is 
set  to  be  from  zero  to  the  tracker  reading  in  axial  direction  padded  in  both  sides 
with  the  maximum  error  of  the  tracker.  Given  two  frame  I\  and  I2 ,  the  amount 
of  lateral  motion  a,  and  compression,  5,  is  found  by  solving  cross-correlation: 


arg  max 

a,b 


E 

x,y£G 


h(x,  y)  ■  I2{x  +  a,  by)  +  h(x  -  a, -by)  ■  I2{x ,  y) 


(4) 


The  RF  data  is  normalized  with  standard  variation  and  assumed  to  have  zero 
mean.  We  employ  two  tricks  which  extensively  increases  the  speed  of  search. 
First,  we  do  not  match  the  entire  image  to  solve  for  these  parameters.  Instead, 
only  pixels  on  a  grid,  G,  are  used  as  described  by  Equation  (4).  The  two  terms  of 
Equation  (4)  ensures  that  the  search  remains  reciprocal,  which  means  switching 
the  images  only  affects  the  sign  of  a  and  b.  Second,  a  is  recovered  by  matching 
only  the  top  part  of  the  two  images  while  b  is  fixed  to  one.  The  reason  is  that 
the  displacement  due  to  compression  is  minimal  in  that  region. 

Having  the  global  motions,  the  cost  function  is  modified  to  penalize  very  low 
compressions: 
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C2(7))  =KX-DX2  +  Ky-^U  topt|3  +Kz-D2z,  (5) 

Dy  +  C 

where  Dx  and  Dy  are  the  global  motions  from  Equation  (4)  converted  to  mm.  c 
is  a  small  number  that  limits  the  cost  of  zero  compression.  Finally,  the  pairs  with 
the  lowest  cost  are  selected  until  a  maximum  number  of  frame  pairs  is  reached 
or  the  minimum  cost  grows  higher  than  the  average  cost. 

The  last  step  involves  computing  the  strain  for  all  the  selected  frame  pairs. 
We  have  implemented  normalized  cross-correlation  (NCC)  [10]  to  recover  the 
displacements  and  least  squares  estimation  to  calculate  the  strain.  Before  calcu¬ 
lating  strain,  the  global  lateral  motion  and  compression  from  the  previous  step 
are  compensated  in  one  image  using  cubic  interpolation.  This  is  known  to  reduce 
the  error  of  strain  estimation  [11].  The  final  strain  image,  Sfinai  is  the  weighted 
average  of  all  the  strains: 


S 


final 


TZi  Wj  -  Sj 

Em 

i=iwi 


Wj  = 


1  Pi 

0, 


Pi  >  0.7 

otherwise 


(6) 


where  pi  is  the  correlation  coefficient  for  the  ith  pair  after  applying  the  dis¬ 
placements,  and  m  is  the  number  of  pairs.  Fusing  the  strains  in  this  fashion  is 
acceptable  since  the  algorithm  only  allows  for  compressions  that  are  close  to  a 
predetermined  amount  optimal  for  strain  estimation. 


3  Experiments  and  Results 

We  acquired  ultrasound  data  using  a  SONOLINE  Antares™  ultrasound  system 
(Siemens  Medical  Solutions  USA,  Inc.)  with  a  high-frequency  ultrasound  trans¬ 
ducer  (VF10-5)  at  center  frequency  of  6-8  MHz.  We  accessed  RF  through  the 
Axius  Direct™  Ultrasound  Research  Interface  provided  by  Siemens.  Our  custom 
data  acquisition  program  was  connected  to  this  interface  to  send  the  command 
for  capturing  RF  data.  At  the  same  time,  the  program  collected  tracking  infor¬ 
mation  from  either  a  “Polaris”  optical  tracker  (Northern  Digital  Inc.,  Waterloo, 
Canada)  with  passive  markers  or  the  “medSAFE”  EM  tracker  (Ascension  Tech. 
Corp.). 

RF  data  and  tracking  information  was  captured  from  a  breast  phantom  con¬ 
taining  a  harder  lesion  (CIRS  elastography  phantom,  Norfolk,  VA)  and  ex-vivo 
pig  liver.  Alginate  was  injected  to  the  liver  to  mark  a  part  of  liver,  and  then, 
that  area  was  ablated.  The  users  were  asked  to  palpate  the  tissue  over  the  hard 
lesion  in  the  breast  phantom  and  the  ablated  lesion  in  the  pig  liver  while  data 
was  being  collected.  Between  100  to  138  RF  frames  were  acquired  with  the  rate 
of  about  30  frames  per  second. 

The  first  set  of  data  was  captured  by  an  experienced  user  from  the  breast 
phantom.  Figure  1(a)  shows  the  translation  components  of  hand  motion  with 
respect  to  the  first  frame.  The  axial  motion  is  dominant  and  there  is  only  a 
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Hand  Motion 


(b)  Strain  image  (TrUE) 


Fig.  1.  (a)  shows  the  translation  of  probe  w.r.t.  the  first  image.  Proper  hand  motion 
is  applied  as  the  axial  compression  is  dominant,  (b)  is  the  output  of  our  proposed 
algorithm. 


gradual  drift  in  the  lateral  and  elevational  directions.  Figure  1(b)  depicts  the 
high-quality  strain  image  resulting  from  the  TrUE  algorithm. 

Applying  a  compression  similar  to  the  one  shown  in  Figure  1(a)  is  a  difficult 
task  for  novice  or  even  intermediate  users.  This  is  especially  the  case  where  axial 
compression  does  not  translate  into  a  simple  up  and  down  motion.  Ultrasound  gel 
creates  a  slippery  surface  that  makes  the  palpation  prone  to  out-of-plane  motion. 
Two  case  are  shown  in  Figure  2,  where  one  is  tracked  with  the  EM  tracker  and 
the  other  one  with  the  optical  tracker.  In  Figure  2(a)  the  hand  motion  contains  a 
large  amount  of  out-of-plane  motion,  whereas,  in  Figure  2(b),  the  user  has  moved 
the  probe  laterally.  In  both  cases,  the  TrUE  algorithm  generates  reliable  results. 
Figures  2  (c)  and  (d)  show  the  contrast-to-noise  ratio  (CNR)  and  signal-to-noise 
ratio  (SNR)  of  the  strain  image.  The  CNR  and  SNR  value  are  computed  from: 


CNR  = 


2(sfa  -  St)2 

vl+°t 


SNR  = 

a 


(7) 


where  s  and  a  denote  the  mean  and  standard  deviation  of  intensities.  The  t  or 
5  subscripts  show  that  the  computation  is  only  for  the  target  or  the  background 
region,  respectively.  The  SNR  and  CNR  for  computing  the  strain  from  consecu¬ 
tive  frames  (the  dashed  curve)  is  compared  to  the  SNR  and  CNR  of  the  strain 
image  from  the  proposed  method  (solid  line).  Using  consecutive  frames  is  the 
standard  method  of  elastography  in  ultrasound  machines.  Almost  in  all  cases 
the  TrUE  algorithm  outperforms  the  consecutive  frames  by  a  large  margin. 

Although  the  SNR  and  CNR  provide  quantitative  measures  to  compare  the 
strain  images,  they  do  not  directly  reflect  the  visual  quality  of  strain.  In  Figure 
3,  we  show  results  of  elastography  using  our  frame  selection  technique  as  well  as 
four  other  strain  images  calculated  from  consecutive  frames.  The  Figure  shows 
the  effects  of  improper  compression  in  consecutive  frames  in  the  strain  image. 
At  the  same  time  our  algorithm  provides  a  single  reliable  strain. 
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Hand  Motion 


Hand  Motion 


(a)  Poor  palpation  with  EM  tacking 


(b)  Poor  palpation  with  optical  tracking 


Fig.  2.  Two  cases  of  improper  motions  are  shown  where  the  hand  motion  suffers  from 
large  lateral  and  elevational  components  evident  in  relative  translations.  The  results 
of  case  1  with  EM  tacker  is  shown  on  the  left  column,  and  the  results  of  case  2  with 
optical  tracker  is  shown  on  the  right  column. 


4  Discussion 

We  presented  a  method  of  ultrasound  elastography  which  is  robust  to  the  quality 
of  the  hand  motion  of  the  user.  Using  the  information  from  an  external  tracker, 
it  automatically  selects  multiple  frame  pairs  with  a  specific  compression  and 
minimal  undesired  motions.  Our  approach  does  not  take  into  account  the  tissue 
motion  from  other  sources  such  as  breathing  or  patient  motion.  However,  these 
types  of  motions  are  not  normally  problematic  since  they  occur  with  a  slower 
pace  compared  to  hand  motion. 

Our  experiments  shows  that  even  when  the  transducer  has  severe  lateral  or 
out-of-plane  motions,  the  algorithm  still  manages  to  produce  good  results.  The 
multi-stage  frame  selection  and  careful  image  fusion  makes  the  TrUE  method 
less  sensitive  to  tacker  accuracy  and  robust  to  strain  estimation  failures. 

We  are  planning  to  use  the  proposed  method  in  a  breast  cancer  study.  For  this 
purpose,  we  will  be  implementing  our  MATLAB  code  in  C.  The  stain  estimation 
which  is  still  the  bottleneck  of  our  approach  will  be  executed  in  GPU  allowing 
for  the  use  of  sophisticated  algorithms. 


P.  Foroughi,  H.  Rivaz,  I.  Fleming,  G.  Hager,  E.  Boctor 


(a)  B-mode  (b)  TrUE  (c)  1V2 


Fig.  3.  Comparison  of  the  strain  from  TrUE  vs.  consecutive  frames  for  ex-vivo  pig  liver. 
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Ultrasound  Elastography  Using  Three  Images 
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Johns  Hopkins  University 


Abstract.  Displacement  1  estimation  is  an  essential  step  for  ultrasound 
elastography  and  numerous  techniques  have  been  proposed  to  improve 
its  quality  using  two  frames  of  ultrasound  RF  data.  This  paper  intro¬ 
duces  a  technique  for  calculating  a  displacement  field  from  three  frames 
of  ultrasound  RF  data.  To  this  end,  we  first  introduce  constraints  on 
variations  of  the  displacement  field  with  time  using  mechanics  of  ma¬ 
terials.  These  constraints  are  then  used  to  generate  a  regularized  cost 
function  that  incorporates  amplitude  similarity  of  three  ultrasound  im¬ 
ages  and  displacement  continuity.  We  optimize  the  cost  function  in  an 
expectation  maximization  (EM)  framework.  Iteratively  reweighted  least 
squares  (IRLS)  is  used  to  minimize  the  effect  of  outliers.  We  show  that, 
compared  to  using  two  images,  the  new  algorithm  reduces  the  noise  of 
the  displacement  estimation.  The  displacement  field  is  used  to  generate 
strain  images  for  quasi-static  elastography.  Phantom  experiments  and 
in-vivo  patient  trials  of  imaging  liver  tumors  and  monitoring  thermal 
ablation  therapy  of  liver  cancer  are  presented  for  validation. 


1  Introduction 

Displacement,  motion  or  time  delay  estimation  in  ultrasound  images  is  an  essen¬ 
tial  step  in  numerous  medical  imaging  tasks  including  the  rapidly  growing  field 
of  imaging  the  mechanical  properties  of  tissue  [1].  In  this  work,  we  perform  dis¬ 
placement  estimation  for  quasi-static  ultrasound  elastography  [1] ,  which  involves 
deforming  the  tissue  slowly  with  an  external  mechanical  force,  imaging  the  tis¬ 
sue  during  the  deformation,  and  performing  displacement  estimation  using  the 
images.  More  specifically,  we  focus  on  real-time  freehand  palpation  elastography 
[2-7]  where  the  external  force  is  applied  by  simply  pressing  the  ultrasound  probe 
against  the  tissue.  Ease  of  use,  real-time  performance  and  providing  invaluable 
elasticity  images  for  diagnosis  and  guidance/monitoring  of  surgical  operations 
are  the  key  factors  that  have  led  to  its  successful  commercialization. 

A  typical  ultrasound  frame  rate  is  20-60  fps.  As  a  result,  an  entire  series  of 
ultrasound  images  are  freely  available  during  the  tissue  deformation.  Multiple 
ultrasound  images  have  been  used  before  to  obtain  strain  images  of  highly  com¬ 
pressed  tissue  by  accumulating  the  intermediate  strain  images,  and  to  obtain 
persistently  high  quality  strain  images  by  performing  weighted  averaging  of  the 
strain  images  [8-10].  Accumulating  and  averaging  strain  images  increases  their 
signal  to  noise  ratio  (SNR)  and  contrast  to  noise  ratio  (CNR).  However,  these 
techniques  are  susceptible  to  drift,  a  problem  with  any  sequential  tracking  system 
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Fig.  1.  Left:  in-vivo  images  of  liver.  First  and  second  (Si  and  S2  from  left)  are  two 
strain  fields  calculated  from  I\  and  /2,  and  from  I2  and  I3  respectively.  Si  &  S2  look 
“similar” .  Third  image  is  Si  —  77S2  for  77  =  1.1.  The  strain  range  in  the  first  two  images 
is  0  to  0.6%,  and  in  the  third  image  is  ±0.3%.  Right  shows  the  ElastMI  algorithm. 


[11].  In  addition,  these  techniques  do  not  exploit  additional  images  to  improve 
displacement  estimation ,  which  has  many  applications  besides  strain  estimation. 
Time  series  of  ultrasound  data  has  also  been  used  to  characterize  tissue  [12]  and 
improve  elasticity  reconstruction  [13]  and  viscoelastic  parameters  [14,  15]. 

Figure  1  shows  two  consecutive  strain  images  calculated  from  three  ultrasound 
images  using  the  2D  analytic  minimization  (AM)  method  [16].  Our  motivation  is 
to  utilize  the  similarity  of  these  two  images  to  calculate  a  low  variance  displace¬ 
ment  field  from  three  images.  The  contributions  of  this  work  are:  (1)  introducing 
constraints  on  variation  of  the  motion  fields  based  on  similarities  of  strain  images 
through  time ;  (2)  proposing  an  EM  algorithm  to  solve  for  motion  fields  using 
three  images,  and  (3)  reporting  clinical  tstudies  of  ablation  guidance/monitoring, 
with  data  collection  corresponding  to  before,  during  and  after  ablation. 

The  rest  of  this  paper  is  summarized  as  follows.  We  first  introduce  the  Elas- 
tography  using  Multiple  Images  (ElastMI)  algorithm  for  tissue  displacement  es¬ 
timation,  which  minimizes  a  cost  function  that  incorporates  data  obtained  from 
three  images  and  exploits  mechanical  constraints.  The  estimated  low  variance 
displacement  field  can  be  used  in  numerous  applications  in  imaging  mechanical 
properties  of  tissue;  we  use  it  for  generating  strain  images  by  calculating  its  spa¬ 
tial  derivative.  We  use  phantom  and  in-vivo  clinical  studies  to  compare  ElastMI 
versus  the  recently  developed  elastography  technique  of  2D  AM  (code  available 
online  at  www.cs.jhu.edu/~rivaz)  [16]. 


2  ElastMI:  Elastography  Using  Multiple  Images 

We  have  a  set  of  p  =  3  images  /&,  k  =  1  •  •  •  3,  each  of  size  m  x  n.  Let  the  2D 
displacement  field  dfc  =  (ak,lk)  denote  the  displacement  between  A  and  ii, 
where  a  refers  to  the  axial  (i.e.  in  the  direction  of  the  ultrasound  beam)  and  l  to 
the  lateral  (i.e.  perpendicular  to  the  beam  and  in  the  imaging  plane)  directions. 
By  the  choice  of  reference  d1  =  0.  Note  that  we  set  I\  as  the  reference  image 
to  simplify  the  notation.  However,  in  our  implementation  we  always  take  the 
middle  image  (i.e.  I 2)  as  the  reference.  Our  goal  is  to  calculate  a  high  quality  d2 
by  utilizing  all  three  images  in  a  group-wise  approach. 
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It  is  well  known  that  many  tissue  types  display  linear  strain-stress  relation  in 
the  0  to  5%  range  (see  [1]  for  example).  In  a  freehand  palpation  elastography 
setup  with  ultrasound  acquisition  rate  of  20  fps  or  more,  taking  three  consecutive 
images  as  7i,  72,  I3  corresponds  to  strain  values  of  less  than  1%  and  therefore 
the  linearity  assumption  is  valid.  Using  this  property  and  some  simplifying  as¬ 
sumptions,  it  can  be  shown  that  the  ratio  of  the  strain  and  displacement  fields 
in  different  times  is  a  constant  value,  i.e.  strain  images  are  similar  up  to  a  scale 
as  in  Figure  1.  We  denote  the  scale  factor  by  77  =  (r]a,  77/),  and  allow  it  to  slightly 
change  spatially  to  account  for  small  nonlinearities  in  the  tissue.  As  such,  ga 
and  rji  are  themselves  scale  fields  in  the  axial  and  lateral  directions  each  of  size 
mxn.  Using  this  notation  we  have  a3  =  r]a.  *a2  and  l3  =  rji.  */2  where  .*  denotes 
point-wise  multiplication. 

Let  6  contain  all  the  displacement  unknowns  d2  and  d3.  The  MAP  estimate 
of  6  is  obtained  by  maximizing  its  posterior  probability 

Pr(0  |  71,J2,I3)  (xPr(hJ2J3  |  0)  Pr(0)  (1) 

where  we  have  ignored  the  normalization  denominator.  The  data  term  is  cal¬ 
culated  as  Pr(/i,  72,^3  |  0)  =  W7  Pr(/i,  72,  ^3,  77  \  0).  The  summation  over  the 
latent  variable  rj  makes  the  optimization  problem  intractable.  We  therefore  use 
Expectation  Maximization  (EM)  to  make  the  problem  tractable  as  following. 

1.  Initialize:  find  an  estimate  for  0  by  applying  the  2D  AM  method  [16]  to  two 
pairs  of  images  (I1J2)  and  (I1J3)  independently. 

2.  E-step:  find  an  estimate  for  77  using  6  (details  below). 

3.  M-step:  update  6  with  the  current  estimate  of  77  (details  below). 

4.  Iterate  between  2  and  3  until  convergence. 

The  algorithm  is  shown  in  Figure  1  right.  Note  that  unlike  the  traditional  EM 
which  maximizes  Pr(/i,  72,73  |  0),  we  maximize  the  posterior  probability  of  0 
(Equation  1).  Steps  2  and  3  are  elaborated  below. 

Calculating  77  from  6  Using  Least  Squares.  At  each  sample  (i,  j)  in  the 

displacement  field  d2j,  i  =  1  •  •  •  m,  j  =  1-  •  -  n  take  a  window  of  size  mw  x  nw 
centered  at  the  sample  (mw  and  nw  are  in  the  axial  and  lateral  directions  re¬ 
spectively  and  both  are  odd  numbers).  Stack  the  axial  and  lateral  components  of 
d2  •  that  are  in  the  window  in  two  vectors  a?  -  and  l2  ■ ,  each  of  length  mw  x  nw . 
Similarly,  generate  a3^  and  l3^  using  d3.  Note  that  since  both  displacement 
fields  d2j  and  d3j  are  calculated  with  respect  to  samples  on  7i,  they  corre¬ 
spond  to  the  same  sample  (i,  j).  We  first  calculate  the  axial  component  r]^j\a 
(77 =  (V(i,j),airl(i,j),i))-  Discarding  the  spatial  information  in  a2 ^  and  a3^-, 
we  can  average  the  two  vectors  into  two  scalers  a2  •  and  a3  -  and  simply  cal¬ 
culate  77 =  a3j/ a?j.  However,  a  more  elegant  way  which  also  takes  into 
account  the  spatial  information  is  by  calculating  the  least  squares  solution  to 
the  following  over-determined  problem  (superscript  T  denotes  transpose). 

a2Ta3  • 

=  a lj  giving  =  JrT  > 

ai,j  ai,j 


(2) 
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which  is  what  we  use  in  our  implementation.  To  calculate  the  ratio  of  the  lateral 
displacement  fields  we  take  into  account  possible  lateral  slip  of  the  probe, 

which  results  in  a  rigid-body-motion.  The  rigid-body-motion  can  be  simply  cal¬ 
culated  by  averaging  the  lateral  displacement  in  d2  ■  and  d3  ^  in  the  entire  image, 
and  calculating  the  difference  between  these  two  average  lateral  displacements. 
The  lateral  scaling  factor  can  be  calculated  using  an  equation  similar  to  2 

where  the  axial  displacement  a ij  is  replaced  with  the  lateral  displacements  1  ij. 
However,  we  use  the  following  approach  which  results  in  a  better  estimate  for 
The  lateral  strain  ei  is  simply  uea  where  v  is  an  unknown  Poisson’s  ratio. 
Since  v  has  a  small  dynamic  range  in  soft  tissue  and  since  the  difference  between 
the  two  displacement  maps  d2  and  d3  is  small,  we  can  assume  that  v  does  not 
vary  from  d2  to  d3.  Therefore,  =  V{i,j),a-  This  gives  better  estimate  for 

since  axial  displacement  estimation  is  more  accurate  [16]. 

Calculating  0  by  Maximizing  Its  Posterior  Probability.  To  analytically 
solve  the  MAP  estimate  of  0,  we  assume  that  the  data  is  independent  and 
that  the  noise  model  is  Gaussian.  Although  not  completely  held  in  real  images, 
these  assumptions  are  also  the  foundation  behind  sum  of  square  difference  and 
correlation  based  elastography  methods,  which  have  been  extensively  shown  to 
produce  reliable  results.  With  these  assumptions,  the  robust  MAP  estimate  for 
6  can  be  obtained  by  minimizing  the  following  cost  function 

m  2 

C(0)  =  E™12>*  (AXi)  -  (Xi  +  d?)  -  5dtT V/2(xj  +  d?)j  + 

i=  1 

rn  2 

E  Wl3’i  (Jl(xi)  “  J3(Xi  +  Vi, ad?)  ~  7?i,a(5dfVJ3(xi  +  J?i,ad|))  + 

i=  1 
rn 

J^idl-dUfAidt-dU)  (3) 

i=  1 

where  d2  is  the  estimate  obtained  using  2D  AM,  £d2  =  d2  —  d2  is  the  update 
in  the  displacement  that  we  are  looking  for,  A  =  diag(a,/3 )  is  a  2  x  2  diagonal 
matrix  with  tunable  regularization  weights  (a,/? )  that  we  adjust  manually  in 
this  work,  and  V  denotes  the  gradient  operator.  Robustness  is  achieved  using 
IRLS  through  weights  Wi2,i  and  ^13,2  which  are  calculated  as  following 

-  f  1  |r-|  <  T 

Wik,i  =  w(/i(xi)  -  4(xj  +  djfe)),for  k  =  2, 3,  and  w(n)  =  <  j_  \\>T  (4) 

l  \ri\  1 7 1  I 

where  T  is  a  tunable  parameter  which  determines  the  residual  level  for  which 
sample  i  can  be  treated  as  outlier.  A  small  T  will  treat  many  samples  as  outliers. 

Setting  the  derivative  of  C  w.r.t.  the  axial  (Saf  =  £d2a)  and  lateral  (£/2  = 
£d2  j)  components  of  £d2  for  i  =  1  •  •  •  m  to  zero  and  stacking  the  2m  unknowns 

in  6d2  =  [ 5a 2  51 2  5a 2  5l\  •  •  •  5a ^  51^] T  and  the  2m  initial  estimates  in  d2  = 

\a\  l\  a2  l\  •  •  •  we  obtain  the  linear  system  of  size  2m: 
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Table  1.  The  SNR  and  CNR  of  the  strain  images  of  Figure  2 


Axial,  2D  AM 

Axial,  ElastMI 

Lateral,  2D  AM 

Lateral,  ElastMI 

SNR 

11.11 

12.64 

6.06 

6.63 

CNR 

8.48 

9.63 

2.96 

3.39 

(Wl2,i  +  Wl3,iVi,a2)I{,a2  («>12  ,i  +  Wl  3,iVi,aVi,l)I[,aI'l,l 

(' Wl2,i  +  ‘Wl3,i‘ni,a‘ni,l)Il,aIl,l  (w12  ,i  +  Wl3,iVi,l2)I[,l 

where  I'2  and  /:'(  are  calculated  respectively  at  fx,  +  d?)  and  at  fx,  +  rji.  *  d?), 
superscript  /  indicates  derivative  and  subscript  a  and  l  determine  whether  the 
derivation  is  in  the  axial  or  lateral  direction,  and  r  is  a  vector  of  length  2 m  with 
elements: 

i  even  :  r*  =  wi2,il'1>a(xi)  h(xi)  -  hfa  +  df)  + 

Wl3,iVi-  *  I[ ,a(xi)  h (x*)  -  h(*i  +  Vi- *  d f) 

i  odd  :  r*  =  w12,il[ti(xi)  h(xi)  -  h (x*  +  df)  + 

Wl3 ,iVi-  *  /l(Xj)  -  J3(Xi  +Vi- *  df)  .  (7) 

The  coefficient  matrix  in  Equation  5  is  pentadiagonal  and  symmetric.  As  such, 
it  can  be  solved  in  8m  operations,  significantly  less  than  (2m)3/3  required  for 
solving  a  full  system.  For  all  the  results  presented  in  this  work,  the  EM  algorithm 
is  iterated  once. 


3  Results  of  Phantom  Experiments  and  Patient  Trials 

RF  data  is  acquired  from  an  Antares  Siemens  system  (Issaquah,  WA)  at  the 
center  frequency  of  6.67  MHz  with  a  VF10-5  linear  array  at  a  sampling  rate 
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width  (mm)  width  (mm)  width  (mm)  width  (mm) 


(a)  2D  AM  axial  (b)  ElastMI  axial  (c)  2D  AM  lateral  (d)  ElastMI  lateral 

Fig.  2.  Strain  images  of  the  CIRS  phantom  with  the  target  and  background  windows 
(for  calculation  of  SNR  and  CNR).  No  Kalman  filter  [16]  is  used  to  ease  the  comparison. 


Table  2.  The  CNR  of  the  strain  images  of  first,  second  and  third  patient  trials  (images 
of  second  patient  are  shown  in  Figure  3).  PI,  P2  and  P3  respectively  correspond  to 
patients  1,  2  and  3.  2 (s&  —  St)2  and  a2  +  a2  indicate  contrast  and  noise  respectively. 


before  ablation 

during  ablation 

after  ablation 

2D  AM 

ElastMI 

2D  AM 

ElastMI 

2D  AM 

ElastMI 

104  x  2(sb  -  St)* 

- 

- 

- 

- 

2.18 

2.22 

PI 

104  x  (a2  +  of) 

- 

- 

- 

- 

0.108 

0.083 

CNR  =  J2(Sb2~%)2 

V  ab+at 

- 

- 

- 

- 

4.49 

5.17 

104  x  2(s6  -  st)z 

0.45 

0.89 

- 

- 

2.08 

2.15 

P2 

104  x  (a2  +  at) 

0.0036 

0.0045 

- 

- 

0.204 

0.142 

CNR  =  ,/2(s"rs"  2)2 

V  ab+c,i 

11.16 

14.05 

- 

- 

3.19 

3.89 

104  x  2(s6  —  st)'2 

0.235 

0.234 

0.0745 

0.1716 

4.85 

4.82 

P3 

104  x  (of  +  a}) 

0.0045 

0.0036 

0.0091 

0.0161 

0.204 

0.171 

CNR  = 

V  ai+<ri 

7.22 

8.01 

2.87 

3.26 

4.87 

5.31 

of  40  MHz.  An  elastography  phantom  (CIRS  elastography  phantom,  Norfolk, 
VA)  is  compressed  axially  in  two  steps  using  a  linear  stage,  and  three  images 
are  acquired.  Resulting  strain  images  are  shown  in  Figure  2.  The  unitless  metric 
signal  to  noise  ratio  (SNR  =  and  contrast  to  noise  ratio  (CNR  =  P~) 

[1]  of  the  ElastMI  algorithm  are  shown  in  Table  1  (The  SNR  is  only  calculated 
for  the  background  window).  Comparing  to  the  2D  AM,  the  ElastMI  algorithm 
improves  the  SNR  by  approximately  14%  and  the  CNR  by  approximately  11%. 
The  high  quality  of  the  lateral  strain  image,  compared  to  state  of  the  art  strain 
imaging  techniques,  is  visually  noticeable. 

In  the  clinical  studies,  RF  data  was  acquired  from  ablation  therapy  of  three 
patients  with  liver  cancer  using  the  Siemens  Antares  ultrasound  machine  in 
the  following  way:  for  the  first  patient  only  after  ablation,  for  the  second  patient 
before  and  after  ablation,  and  for  the  third  patient  before,  during  and  after  abla¬ 
tion.  The  ablation  was  administered  using  the  RITA  Model  1500  XRF  generator 
(Rita  Medical  Systems,  Fremont,  CA).  Tissue  was  simply  compressed  freehand 
at  a  frequency  of  approximately  1  compression  per  2  sec  with  the  ultrasound 
probe  without  any  attachment,  and  the  strain  images  are  generated  offline. 
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Fig.  3.  Axial  strain  images  of  the  second  in-vivo  patient  trial  corresponding  to  before 
(top  row)  and  after  (bottom  row)  ablation.  The  first,  second  and  third  columns  are 
respectively  B-mode,  2D  AM  strain  and  ElastMI  strain  images.  The  cancer  tumor  in 
the  top  row,  and  the  ablated  lesion  in  the  bottom  row  are  delineated.  The  CNR  between 
the  target  and  background  (marked  by  t  &  b)  windows  are  given  in  Table  2. 


Results  of  the  second  patient  trial  are  shown  in  Figure  3.  Considering  the 
numerous  sources  of  noise  in  the  clinical  data,  the  high  contrast  of  the  tumor 
(top  row)  and  the  ablated  lesion  (bottom  row)  in  the  strain  images  make  ElastMI 
a  promising  tool  for  both  finding  the  tumor  and  monitoring  the  ablation.  It 
should  be  noted  that  elastographic  analysis  of  the  ablated  lesion  is  known  to  be 
challenging  due  to  high  temperatures  which  significantly  degrade  the  quality  of 
ultrasound  data  (mainly  because  of  the  air  bubbles).  Table  2  summarizes  the 
CNR,  as  well  as  noise  and  contrast  values,  in  the  patient  trials  obtained  using 
2D  AM  and  ElastMI  methods.  In  the  six  cases  presented  in  this  table  (two  before 
ablation,  one  during  ablation  and  three  after  ablation),  the  average  increase  in 
the  CNR  achieved  using  ElastMI  compared  to  2D  AM  is  17%. 


4  Conclusions 

In  this  work,  we  proposed  to  utilize  three  ultrasound  images  to  calculate  high 
quality  displacement  fields.  We  neglected  the  dynamics  of  tissue  motion  and 
assumed  a  static  model  for  tissue  mechanics,  which  is  valid  in  the  quasi-static 
elastography.  Using  this  model  and  assuming  tissue  linearity,  which  holds  in 
the  low  strain  rates  of  the  freehand  elastography,  we  introduced  constraints  on 
the  variations  of  the  strain  field  with  time.  We  then  proposed  ElastMI,  an  EM 
algorithm  that  exploits  these  constraints  for  estimating  displacement  fields  using 
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three  images.  The  algorithm  involves  solving  sparse  linear  systems,  and  therefore 
runs  in  real-time.  The  low  variance  motion  field  that  we  compute  by  exploiting 
this  new  prior  can  be  used  in  numerous  applications  in  ultrasound  imaging;  we 
used  it  here  to  generate  strain  images. 
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